1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
include geoplot_utils.inp
# public
function bundle geoplot_describe_json (const bundle jb, int verbose[1])
# as per RFC 7946, "The GeoJSON Format"
bundle ret = null
matrix bbox = get_bbox(jb)
if rows(bbox) == 0
printf "No bounding box\n"
else
printf "%d-dimensional bounding box\n", rows(bbox)
endif
ndim = rows(bbox)
ret.bbox = bbox
bundles feat = jb.features
nf = nelem(feat)
printf "%d features\n", nf
ret.nf = nf
bundle ids = null
loop i = 1 .. nf
props = feat[i].properties
if i == 1
ids = create_ids(props, nf)
endif
fill_ids(props, &ids, i)
if verbose
describe_feature(feat[i], i, 0)
endif
endloop
ret.ids = ids
return ret
end function
# public
function void geoplot_translate_feature(bundle *b, int f,
matrix shift,
matrix center[null],
matrix scale[null])
if exists(center)
center = vec(center)
else
center = {0;0}
endif
if exists(scale)
scale = vec(scale)
else
scale = {1;1}
endif
shift = vec(shift)
scalar single = b.features[f].geometry.type == "Polygon"
arrays mod = b.features[f].geometry.coordinates
loop i=1..nelem(mod)
loop j=1..nelem(mod[i])
if single
ctmp = scale .* (mod[i][j] - center)
mod[i][j] = ctmp + center + shift
else
loop k=1..nelem(mod[i][j])
ctmp = scale .* (mod[i][j][k] - center)
mod[i][j][k] = ctmp + center + shift
endloop
endif
endloop
endloop
b.features[f].geometry.coordinates = mod
end function
function void geoplot_set_properties (bundle *b, list L)
strings keys = varnames(L)
loop i=1..nelem(b.features)
bundle bpi = null
loop foreach j L
bpi[keys[j]] = L.$j[i]
endloop
b.features[i].properties = bpi
endloop
end function
function matrix geoplot_seek_feature(const bundle b,
string name,
bool do_plot[1])
bundles feat = b.features
nf = nelem(feat)
matrix found = {}
scalar nfound = 0
strings where = array(0)
loop i = 1 .. nf
bundle prop_i = feat[i].properties
strings keys = sort(getkeys(prop_i))
loop j = 1 .. nelem(keys)
type = inbundle(prop_i, keys[j])
if type == 4
if instring(tolower(prop_i[keys[j]]), tolower(name))
found = found ~ i
nfound++
where = where + defarray(keys[j])
break
endif
endif
endloop
endloop
if nfound == 0
printf "\"%s\" not found\n", name
elif nfound == 1
f = feat[found[1]]
extents = describe_feature(f, 0, 0)
string key = where[1]
string cont = f.properties[key]
if do_plot
plot_feature(b, found[1], cont, extents)
endif
else
loop i = 1 .. nfound
k = found[i]
string key = where[i]
string cont = feat[k].properties[key]
printf "Feature %4d: \"%s\" = %s\n", found[i], key, cont
endloop
endif
return found
end function
function void geoplot_simplify(bundle *b,
scalar preserve[0.1:1:0.75])
feat = b.features
set stopwatch
loop i = 1 .. nelem(feat)
feat[i] = simplify_feature(feat[i], preserve)
endloop
printf "Simplification took %g seconds\n", $stopwatch
flush
b.features = feat
end function
|
