I was hoping to create an equidistant grid over Manhattan map (say 200m by 200m) projection using latitudes and longitudes not degrees. I was using basemap but couldn't figure out a way to proceed with the task. This the projection code.
m = Basemap(projection='mill',
llcrnrlat= 40.6968,
llcrnrlon= -74.0224,
urcrnrlat= 40.8964,
urcrnrlon= -73.8927,
resolution='h')
What will be the best way to do the above, I also need to store lat,long values of each grid vertex points.
From the basemap documentation:
In order to plot data on a map, the coordinates of the data must be given in map projection coordinates. Calling a Basemap class instance with the arguments lon, lat will convert lon/lat (in degrees) to x/y map projection coordinates (in meters). The inverse transformation is done if the optional keyword inverse is set to True.
There is also an example in the documentation page. To adept this example to your use case, I converted the lower left corners into meters, produced a regular grid with 2000 m spacing (200 m was a bit too dense) and converted the grid back into lon/lat coordinates, which can then be used by drawparallels()
and drawmeridians
.
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
fig, ax = plt.subplots()
lon0 = -74.0224
lat0 = 40.6968
m = Basemap(
projection='mill',
llcrnrlat= 40.6968,
llcrnrlon= -74.0224,
urcrnrlat= 40.8964,
urcrnrlon= -73.8927,
resolution='h',
ax = ax,
)
x0, y0 = m(lon0, lat0)
x = np.arange(20)*2000+x0
y = np.arange(20)*2000+y0
lons, lats = m(x,y,inverse=True)
m.drawcoastlines()
m.drawcountries()
m.drawstates()
m.drawmeridians(lons)
m.drawparallels(lats)
plt.show()
The result looks like this:
Hope this helps.