1#!/usr/bin/env python3
2###############################################################################
3# $Id: mkgraticule.py ffccab1ee20b5151e9bd45f1a2c46245a74f1f56 2021-04-23 14:05:42 +0300 Idan Miara $
4#
5# Project:  OGR Python samples
6# Purpose:  Produce a graticule (grid) dataset.
7# Author:   Frank Warmerdam, warmerdam@pobox.com
8#
9###############################################################################
10# Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
11#
12# Permission is hereby granted, free of charge, to any person obtaining a
13# copy of this software and associated documentation files (the "Software"),
14# to deal in the Software without restriction, including without limitation
15# the rights to use, copy, modify, merge, publish, distribute, sublicense,
16# and/or sell copies of the Software, and to permit persons to whom the
17# Software is furnished to do so, subject to the following conditions:
18#
19# The above copyright notice and this permission notice shall be included
20# in all copies or substantial portions of the Software.
21#
22# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28# DEALINGS IN THE SOFTWARE.
29###############################################################################
30
31import sys
32
33from osgeo import ogr
34from osgeo import osr
35
36
37#############################################################################
38def float_range(*args):
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#############################################################################
61def Usage():
62    print('Usage: mkgraticule [-connected] [-s stepsize] [-substep substepsize]')
63    print('         [-t_srs srs] [-range xmin ymin xmax ymax] outfile')
64    print('')
65    return 1
66
67
68def main(argv):
69    # Argument processing.
70    t_srs = None
71    stepsize = 5.0
72    substepsize = 5.0
73    connected = 0
74    outfile = None
75
76    xmin = -180
77    xmax = 180
78    ymin = -90
79    ymax = 90
80
81    i = 1
82    while i < len(argv):
83        if argv[i] == '-connected':
84            connected = 1
85        elif argv[i] == '-t_srs':
86            i = i + 1
87            t_srs = argv[i]
88        elif argv[i] == '-s':
89            i = i + 1
90            stepsize = float(argv[i])
91        elif argv[i] == '-substep':
92            i = i + 1
93            substepsize = float(argv[i])
94        elif argv[i] == '-range':
95            xmin = float(argv[i + 1])
96            ymin = float(argv[i + 2])
97            xmax = float(argv[i + 3])
98            ymax = float(argv[i + 4])
99            i = i + 4
100        elif argv[i][0] == '-':
101            return Usage()
102        elif outfile is None:
103            outfile = argv[i]
104        else:
105            return Usage()
106
107        i = i + 1
108
109    if outfile is None:
110        outfile = "graticule.shp"
111
112
113    if substepsize > stepsize:
114        substepsize = stepsize
115
116    # -
117    # Do we have an alternate SRS?
118
119    ct = None
120
121    if t_srs is not None:
122        t_srs_o = osr.SpatialReference()
123        t_srs_o.SetFromUserInput(t_srs)
124
125        s_srs_o = osr.SpatialReference()
126        s_srs_o.SetFromUserInput('WGS84')
127
128        ct = osr.CoordinateTransformation(s_srs_o, t_srs_o)
129    else:
130        t_srs_o = osr.SpatialReference()
131        t_srs_o.SetFromUserInput('WGS84')
132
133    # -
134    # Create graticule file.
135
136    drv = ogr.GetDriverByName('ESRI Shapefile')
137
138    try:
139        drv.DeleteDataSource(outfile)
140    except:
141        pass
142
143    ds = drv.CreateDataSource(outfile)
144    layer = ds.CreateLayer('out', geom_type=ogr.wkbLineString,
145                           srs=t_srs_o)
146
147    #########################################################################
148    # Not connected case.  Produce individual segments are these are going to
149    # be more resilient in the face of reprojection errors.
150
151    if not connected:
152        #########################################################################
153        # Generate lines of latitude.
154
155        feat = ogr.Feature(feature_def=layer.GetLayerDefn())
156        geom = ogr.Geometry(type=ogr.wkbLineString)
157
158        for lat in float_range(ymin, ymax + stepsize / 2, stepsize):
159            for long_ in float_range(xmin, xmax - substepsize / 2, substepsize):
160
161                geom.SetPoint(0, long_, lat)
162                geom.SetPoint(1, long_ + substepsize, lat)
163
164                err = 0
165                if ct is not None:
166                    err = geom.Transform(ct)
167
168                if err == 0:
169                    feat.SetGeometry(geom)
170                    layer.CreateFeature(feat)
171
172        #########################################################################
173        # Generate lines of longitude
174
175        for long_ in float_range(xmin, xmax + stepsize / 2, stepsize):
176            for lat in float_range(ymin, ymax - substepsize / 2, substepsize):
177                geom.SetPoint(0, long_, lat)
178                geom.SetPoint(1, long_, lat + substepsize)
179
180                err = 0
181                if ct is not None:
182                    err = geom.Transform(ct)
183
184                if err == 0:
185                    feat.SetGeometry(geom)
186                    layer.CreateFeature(feat)
187
188
189    #########################################################################
190    # Connected case - produce one polyline for each complete line of latitude
191    # or longitude.
192
193    if connected:
194        #########################################################################
195        # Generate lines of latitude.
196
197        feat = ogr.Feature(feature_def=layer.GetLayerDefn())
198
199        for lat in float_range(ymin, ymax + stepsize / 2, stepsize):
200
201            geom = ogr.Geometry(type=ogr.wkbLineString)
202
203            for long_ in float_range(xmin, xmax + substepsize / 2, substepsize):
204                geom.AddPoint(long_, lat)
205
206            err = 0
207            if ct is not None:
208                err = geom.Transform(ct)
209
210            if err == 0:
211                feat.SetGeometry(geom)
212                layer.CreateFeature(feat)
213
214        #########################################################################
215        # Generate lines of longitude
216
217        for long_ in float_range(xmin, xmax + stepsize / 2, stepsize):
218
219            geom = ogr.Geometry(type=ogr.wkbLineString)
220
221            for lat in float_range(ymin, ymax + substepsize / 2, substepsize):
222                geom.AddPoint(long_, lat)
223
224            err = 0
225            if ct is not None:
226                err = geom.Transform(ct)
227
228            if err == 0:
229                feat.SetGeometry(geom)
230                layer.CreateFeature(feat)
231
232    #############################################################################
233    # Cleanup
234
235    feat = None
236    geom = None
237
238    ds.Destroy()
239    ds = None
240
241    return 0
242
243
244if __name__ == '__main__':
245    sys.exit(main(sys.argv))
246