I am trying to create a fishnet of x by x miles for all USA, same as this basically: http://www.giscourse.com/creating-a-fishnet-grid-using-arcgis-10/
However, I want to do this with python.
My current code works but it doesn't take distance in miles, but rather in degrees, and therefore it doesn't take into consideration how for different lat, lon points the conversion of distance to degrees changes (Calculating longitude length in miles?).
This is my current code:
from pygeoif import geometry
import shapefile
import numpy as np
sq_len_deg = 1
lats = np.arange(26, 50.1, sq_len_deg)
lons = np.arange(-125, -69.9, sq_len_deg)
a=[]
geometries=[]
for lat in lats:
for lon in lons:
a.append([(lon,lat),(lon,lat-sq_len_deg),(lon-sq_len_deg,lat-sq_len_deg),(lon-sq_len_deg,lat)])
geometries.append('POLYGON ((%s %s, %s %s, %s %s, %s %s))'%(lon,lat,lon,lat-sq_len_deg,lon-sq_len_deg,lat-sq_len_deg,lon-sq_len_deg,lat))
w = shapefile.Writer(shapefile.POLYGON)
w.field('test','N')
for i,poly in enumerate(geometries):
w.poly([geometry.from_wkt(poly).exterior.coords])
print(i)
w.record(i)
w.save('officialpysh')
I can see this going into two directions, some function already exist and I can use this, or a suggestion on how to modify this code.