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