Package osgeo :: Package utils :: Module mkgraticule
[hide private]
[frames] | no frames]

Source Code for Module osgeo.utils.mkgraticule

  1  #!/usr/bin/env python3 
  2  ############################################################################### 
  3  # $Id$ 
  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   
 31  import sys 
 32   
 33  from osgeo import ogr 
 34  from osgeo import osr 
 35   
 36   
 37  ############################################################################# 
38 -def 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 #############################################################################
61 -def Usage():
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
69 -def main(argv):
70 # Argument processing. 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 # Do we have an alternate SRS? 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 # Create graticule file. 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 # Not connected case. Produce individual segments are these are going to 150 # be more resilient in the face of reprojection errors. 151 152 if not connected: 153 ######################################################################### 154 # Generate lines of latitude. 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 # Generate lines of longitude 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 # Connected case - produce one polyline for each complete line of latitude 192 # or longitude. 193 194 if connected: 195 ######################################################################### 196 # Generate lines of latitude. 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 # Generate lines of longitude 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 # Cleanup 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