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 import sys
32
33 from osgeo import ogr
34 from osgeo import osr
35
36
37
39 start = 0.0
40 step = 1.0
41 if len(args) == 1:
42 (stop,) = args
43 elif len(args) == 2:
44 (start, stop) = args
45 elif len(args) == 3:
46 (start, stop, step) = args
47 else:
48 raise TypeError("float_range needs 1-3 float arguments")
49
50 the_range = []
51 steps = (stop - start) / step
52 if steps != int(steps):
53 steps = steps + 1.0
54 for i in range(int(steps)):
55 the_range.append(i * step + start)
56
57 return the_range
58
59
60
62 print('Usage: mkgraticule [-connected] [-s stepsize] [-substep substepsize]')
63 print(' [-t_srs srs] [-range xmin ymin xmax ymax] outfile')
64 print('')
65 sys.exit(1)
66
67
68
70
71 t_srs = None
72 stepsize = 5.0
73 substepsize = 5.0
74 connected = 0
75 outfile = None
76
77 xmin = -180
78 xmax = 180
79 ymin = -90
80 ymax = 90
81
82 i = 1
83 while i < len(argv):
84 if argv[i] == '-connected':
85 connected = 1
86 elif argv[i] == '-t_srs':
87 i = i + 1
88 t_srs = argv[i]
89 elif argv[i] == '-s':
90 i = i + 1
91 stepsize = float(argv[i])
92 elif argv[i] == '-substep':
93 i = i + 1
94 substepsize = float(argv[i])
95 elif argv[i] == '-range':
96 xmin = float(argv[i + 1])
97 ymin = float(argv[i + 2])
98 xmax = float(argv[i + 3])
99 ymax = float(argv[i + 4])
100 i = i + 4
101 elif argv[i][0] == '-':
102 Usage()
103 elif outfile is None:
104 outfile = argv[i]
105 else:
106 Usage()
107
108 i = i + 1
109
110 if outfile is None:
111 outfile = "graticule.shp"
112
113
114 if substepsize > stepsize:
115 substepsize = stepsize
116
117
118
119
120 ct = None
121
122 if t_srs is not None:
123 t_srs_o = osr.SpatialReference()
124 t_srs_o.SetFromUserInput(t_srs)
125
126 s_srs_o = osr.SpatialReference()
127 s_srs_o.SetFromUserInput('WGS84')
128
129 ct = osr.CoordinateTransformation(s_srs_o, t_srs_o)
130 else:
131 t_srs_o = osr.SpatialReference()
132 t_srs_o.SetFromUserInput('WGS84')
133
134
135
136
137 drv = ogr.GetDriverByName('ESRI Shapefile')
138
139 try:
140 drv.DeleteDataSource(outfile)
141 except:
142 pass
143
144 ds = drv.CreateDataSource(outfile)
145 layer = ds.CreateLayer('out', geom_type=ogr.wkbLineString,
146 srs=t_srs_o)
147
148
149
150
151
152 if not connected:
153
154
155
156 feat = ogr.Feature(feature_def=layer.GetLayerDefn())
157 geom = ogr.Geometry(type=ogr.wkbLineString)
158
159 for lat in float_range(ymin, ymax + stepsize / 2, stepsize):
160 for long_ in float_range(xmin, xmax - substepsize / 2, substepsize):
161
162 geom.SetPoint(0, long_, lat)
163 geom.SetPoint(1, long_ + substepsize, lat)
164
165 err = 0
166 if ct is not None:
167 err = geom.Transform(ct)
168
169 if err == 0:
170 feat.SetGeometry(geom)
171 layer.CreateFeature(feat)
172
173
174
175
176 for long_ in float_range(xmin, xmax + stepsize / 2, stepsize):
177 for lat in float_range(ymin, ymax - substepsize / 2, substepsize):
178 geom.SetPoint(0, long_, lat)
179 geom.SetPoint(1, long_, lat + substepsize)
180
181 err = 0
182 if ct is not None:
183 err = geom.Transform(ct)
184
185 if err == 0:
186 feat.SetGeometry(geom)
187 layer.CreateFeature(feat)
188
189
190
191
192
193
194 if connected:
195
196
197
198 feat = ogr.Feature(feature_def=layer.GetLayerDefn())
199
200 for lat in float_range(ymin, ymax + stepsize / 2, stepsize):
201
202 geom = ogr.Geometry(type=ogr.wkbLineString)
203
204 for long_ in float_range(xmin, xmax + substepsize / 2, substepsize):
205 geom.AddPoint(long_, lat)
206
207 err = 0
208 if ct is not None:
209 err = geom.Transform(ct)
210
211 if err == 0:
212 feat.SetGeometry(geom)
213 layer.CreateFeature(feat)
214
215
216
217
218 for long_ in float_range(xmin, xmax + stepsize / 2, stepsize):
219
220 geom = ogr.Geometry(type=ogr.wkbLineString)
221
222 for lat in float_range(ymin, ymax + substepsize / 2, substepsize):
223 geom.AddPoint(long_, lat)
224
225 err = 0
226 if ct is not None:
227 err = geom.Transform(ct)
228
229 if err == 0:
230 feat.SetGeometry(geom)
231 layer.CreateFeature(feat)
232
233
234
235
236 feat = None
237 geom = None
238
239 ds.Destroy()
240 ds = None
241
242
243 if __name__ == '__main__':
244 sys.exit(main(sys.argv))
245