Map plot grid issues

Question

I want to plot a map using Basemap and I can't find how to plot a UTM grid on the map.

I've seen how to plot the grid using long/lat but not in UTM. In Basemap y use epsg=5520 which is UTM 31N.

m = Basemap(epsg=5520, llcrnrlat=52, llcrnrlon=5,urcrnrlat=53, urcrnrlon=6, resolution='l')

m.arcgisimage(server='http://server.arcgisonline.com/arcgis',
              service='World_Imagery', xpixels=3500)
m.drawparallels(np.arange(52, 53, 0.05), labels=[1, 0, 0, 0])
m.drawmeridians(np.arange(5, 6, 0.05), labels=[0, 0, 0, 1])

Any thoughts about how to implement a UTM grid?

score 2 · Accepted Answer · answered Feb 10 '19 at 14:25

With Basemap, plotting UTM grid lines or ticks on UTM map projection is not easy because Basemap's data coordinates (conversion from long-lat) are deviated from real UTM values. So, to get appropriate (x,y) from (long, lat), I use pyproj package. In the provided code, command plot() is used to plot all the grid ticks. And annotate() is used to plot the grid labels outside the map area. Values of grid labels need to multiply with 10000 to get metres units.

Here is the working code:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import pyproj

# need pyproj package for coordinate tranformation
pp = pyproj.Proj(init='epsg:5520')

# use map's corners (long,lat) to get grid coordinates (x,y)
corners = [[5,52], [5,53], [6.2, 53], [5.95, 51.5]]
for ea in corners:
    x,y = pp(ea[0],ea[1])  #(long,lat) to (x,y)
    lon,lat = pp(x, y, inverse=True)
    print(x, y, "%4.1f"%(lon), "%4.1f"%(lat))

# the output of print() above, give extents in grid coordinates
#x range: 1630000, 1720000 m -> 1630, 1720 km
#y range: 5710000, 5900000 m -> 5710, 5900 km

low_x, hi_x = 1630000, 1720000
low_y, hi_y = 5710000, 5900000
grid_sp = 10000   # 10km grid spacing

# we will plot grid ticks '+' at 10km spacing
# .. inside the plotting area
lon_lat = []   # for positions of grid ticks '+'
for ea in np.arange(low_x, hi_x, grid_sp):      # xs
    for eb in np.arange(low_y, hi_y, grid_sp):  # ys
        lon,lat = pp(ea, eb, inverse=True)
        #print(ea, eb, lon, lat)
        # lon, lat is good for plotting on basemap
        lon_lat.append([lon,lat])

# for annotation above top edge, every 10km
yt = 5870000   # y at top edge of map
xs_top = []    # for labels' positions of x grid
for xi in np.arange(low_x, hi_x, grid_sp):   # xs
    lon,lat = pp(xi, yt, inverse=True)
    #print(ea, eb, lon, lat)
    xs_top.append([lon,lat])

# make anno text for every 10 km along map's top edge
anno_top = map(str, list(range(low_x/grid_sp, hi_x/grid_sp)))

# for annotation to the right, every 10km
xr = 1700000   # x at the right edge of map
ys_rt = []     # for labels' positions of y grid
for yi in np.arange(low_y, hi_y, grid_sp):   # ys
    lon,lat = pp(xr, yi, inverse=True)
    #print(xr, yi, lon, lat)
    ys_rt.append([lon,lat])

# make anno text for every 10 km along map's right edge
anno_rt = map(str, list(range(low_y/grid_sp, hi_y/grid_sp))) 

# prep fig/axes for Basemap plot
fig, ax = plt.subplots(figsize=(10, 12))

m = Basemap(epsg=5520, llcrnrlat=52, llcrnrlon=5, urcrnrlat=53, urcrnrlon=6, resolution='i')

# option to plot imagery, need internet connection
if True:
    server = 'http://server.arcgisonline.com/arcgis'
    m.arcgisimage(server=server, service='World_Imagery', xpixels=1500)

# plot grid ticks '+' inside map area
m.plot(np.array(lon_lat)[:,0], np.array(lon_lat)[:,1], 'w+', latlon=True, zorder=10)

# option to plot grid labels on top/right edges
if True:
    # grid labels on top edge
    for id,ea in enumerate(xs_top):
        if ea[0]>5.0 and ea[0]<6.0:
            ax.annotate(anno_top[id], \
                        m(ea[0], ea[1]), \
                        xytext=[-8,50], \
                        textcoords='offset points', \
                        color='b')
            pass
        pass

    # grid labels on right edge
    for id,ea in enumerate(ys_rt):
        if ea[1]>52.0 and ea[1]<53.0:
            ax.annotate(anno_rt[id], \
                        m(ea[0], ea[1]), \
                        xytext=[10,-5], \
                        textcoords='offset points', \
                        color='b')
            pass
        pass

#m.drawcoastlines(linewidth=0.25)
#m.fillcontinents(color='lightgray')

m.drawparallels(np.arange(52, 53.1, 0.1), labels=[1, 0, 0, 0])
m.drawmeridians(np.arange(5, 6, 0.1), labels=[0, 0, 0, 1])
plt.show()

The resulting plot:

Map plot grid issues

1 Answers1