Проблемы с сеткой карты

Question

Я хочу построить карту, используя Базовую карту, и не могу найти способ построения сетки UTM на карте.

Я видел, как построить график сетки, используя long / lat, но не в UTM.В базовой карте y используйте epsg = 5520, что является 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])

Есть мысли о том, как реализовать сетку UTM?

Answer 1

С помощью базовой карты построение линий или отметок сетки UTM на проекции карты UTM не является простым, поскольку координаты данных базовой карты (преобразование из long-lat) отклоняются от реальных значений UTM.Итак, чтобы получить соответствующие (x, y) из (long, lat), я использую пакет pyproj.В предоставленном коде команда plot() используется для построения всех отметок сетки.И annotate() используется для построения меток сетки за пределами области карты.Значения меток сетки нужно умножить на 10000, чтобы получить metres единиц.

Вот рабочий код:

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()

Полученный график: