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