Following Creating square grid polygon shapefile with Python?, I'd like to create an UTM or MGRS graticule in python. In the mentioned question the given answer gives me a way to generate a grid but I don't know if the SRS or the projection is wrong, but the grid is not aligned with the graticule available on NGA.
Here is what I get :
Red grid is the grid from NGA, it aligns well in regard with UTM zone 31 and distance is good too. Blue grid is the one generated by the script I got from Creating square grid polygon shapefile with Python?, but modified to take geographic bounds, and intended utm zone. There is a slight offset (better than previously)
import os, sys
import ogr
import osr
from math import ceil
def main(outputGridfn,lonmin,lonmax,latmin,latmax, utmZone, gridHeight,gridWidth):
# create the spatial reference, WGS84
targetSRS = osr.SpatialReference()
targetSRS.SetProjCS( "UTM "+ utmZone+ " (WGS84) in northern hemisphere." );
targetSRS.SetWellKnownGeogCS( "WGS84" );
targetSRS.SetUTM( int(utmZone), True );
sourceSRS = osr.SpatialReference()
sourceSRS.ImportFromEPSG(4326)
poCT =osr.CoordinateTransformation ( sourceSRS, targetSRS )
# convert sys.argv to float
lonmin = float(lonmin)
lonmax = float(lonmax)
latmin = float(latmin)
latmax = float(latmax)
xmin,ymin,_ = poCT.TransformPoint(lonmin,latmin)
print (xmin, ymin)
xmax,ymax,_ = poCT.TransformPoint(lonmax,latmax)
gridWidth = float(gridWidth)
gridHeight = float(gridHeight)
# get rows
rows = ceil((ymax-ymin)/gridHeight)
# get columns
cols = ceil((xmax-xmin)/gridWidth)
# start grid cell envelope
ringXleftOrigin = xmin
ringXrightOrigin = xmin + gridWidth
ringYtopOrigin = ymax
ringYbottomOrigin = ymax-gridHeight
# create output file
outDriver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(outputGridfn):
os.remove(outputGridfn)
outDataSource = outDriver.CreateDataSource(outputGridfn)
outLayer = outDataSource.CreateLayer(outputGridfn,targetSRS,geom_type=ogr.wkbPolygon )
featureDefn = outLayer.GetLayerDefn()
# create grid cells
countcols = 0
while countcols < cols:
countcols += 1
# reset envelope for rows
ringYtop = ringYtopOrigin
ringYbottom =ringYbottomOrigin
countrows = 0
while countrows < rows:
countrows += 1
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(ringXleftOrigin, ringYtop)
ring.AddPoint(ringXrightOrigin, ringYtop)
ring.AddPoint(ringXrightOrigin, ringYbottom)
ring.AddPoint(ringXleftOrigin, ringYbottom)
ring.AddPoint(ringXleftOrigin, ringYtop)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
# add new geom to layer
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(poly)
outLayer.CreateFeature(outFeature)
outFeature.Destroy
# new envelope for next poly
ringYtop = ringYtop - gridHeight
ringYbottom = ringYbottom - gridHeight
# new envelope for next poly
ringXleftOrigin = ringXleftOrigin + gridWidth
ringXrightOrigin = ringXrightOrigin + gridWidth
# Close DataSources
outDataSource.Destroy()
if __name__ == "__main__":
#
# example run : $ python grid.py <full-path><output-shapefile-name>.shp xmin xmax ymin ymax gridHeight gridWidth
#
if len( sys.argv ) != 9:
print ("[ ERROR ] you must supply eight arguments: output-shapefile-name.shp lonmin lonmax latmin latmax utmZone gridHeight gridWidth")
sys.exit( 1 )
main( sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7],sys.argv[8] )
I'm quite confused at this time. Should I use directly UTM or MGRS coordinates?
I've seen in Generating shapefile containing UTM grid (graticule)? the option of using Saga GIS, but if possible I'd rather use code to be able to script it.
Edit: improved the computation