Skip to main content
 首页 » 编程设计

python之绘制正方形 Cartopy map

2024年09月07日60dudu

我需要使用 Cartopy 绘制方形 map 。我目前为我的 map 使用以下代码:

plt.figure(figsize = (15, 15)) 
 
img = cimgt.GoogleTiles() 
 
ax = plt.axes(projection = img.crs) 
ax.set_extent((d['longitude'].min() - 0.05, d['longitude'].max() + 0.05, 
               d['latitude'].min() - 0.05, d['latitude'].max() + 0.05)) 
 
ax.add_image(img, 10, interpolation = 'bicubic') 
 
plt.scatter(d['longitude'], d['latitude'], transform = ccrs.PlateCarree(), 
            c = '#E8175D', s = 14) 

这工作正常,除了 map 不是正方形这一事实。相反,它正好适合绘图的 (15, 15) 正方形。

我想在左侧和右侧添加更多 map ,使绘图完全呈正方形而不扭曲。简单地将范围设置为相同的纬度和经度差异并不能完成这项工作,因为纬度和经度覆盖 Google(和大多数其他) map 投影中的不同距离。我还找到了this发布,但据我所知,这里的目的是扭曲 map 。

我希望有人知道如何做到这一点。 Cartopy在这方面似乎不是很直观。

请您参考如下方法:

要获得方形范围,您需要在 map 投影坐标中指定它。这涉及到一些坐标变换。这是您需要的代码片段。

# crs of your choice 
crg = cimgt.StamenTerrain().crs    # or cimgt.GoogleTiles().crs 
 
# set map limits, in degrees 
lonmin, lonmax = -22, -15 
latmin, latmax = 63, 65 
 
# do coordinate transformation 
LL = crg.transform_point(lonmin, latmin, ccrs.Geodetic()) 
UR = crg.transform_point(lonmax, latmax, ccrs.Geodetic()) 
EW = UR[0] - LL[0] 
SN = UR[1] - LL[1] 
 
# get side of the square extent (in map units, usually meters) 
side = max(EW, SN)    # larger value is in effect 
mid_x, mid_y = LL[0]+EW/2.0, LL[1]+SN/2.0  # center location 
 
# the extent preserves the center location 
extent = [mid_x-side/2.0, mid_x+side/2.0, mid_y-side/2.0, mid_y+side/2.0] 
 
# this sets square extent 
# crs=crg signifies that projection coordinates is used in extent 
ax.set_extent(extent, crs=crg) 

希望对您有所帮助。

编辑

这是一个完整的工作代码及其生成的 map 。

import matplotlib.pyplot as plt 
import cartopy.crs as ccrs 
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER 
import cartopy.io.img_tiles as cimgt 
 
def make_map(projection=ccrs.PlateCarree()): 
    fig, ax = plt.subplots(figsize=(10, 10), 
                           subplot_kw=dict(projection=projection)) 
    gl = ax.gridlines(draw_labels=True) 
    gl.xlabels_top = gl.ylabels_right = False 
    gl.xformatter = LONGITUDE_FORMATTER 
    gl.yformatter = LATITUDE_FORMATTER 
    return fig, ax 
 
request = cimgt.StamenTerrain()   # very responsive 
crg = request.crs   #crs of the projection 
fig, ax = make_map(projection = crg) 
 
# specify map extent here 
lonmin, lonmax = -22, -15 
latmin, latmax = 63, 65 
 
LL = crg.transform_point(lonmin, latmin, ccrs.Geodetic()) 
UR = crg.transform_point(lonmax, latmax, ccrs.Geodetic()) 
EW = UR[0] - LL[0] 
SN = UR[1] - LL[1] 
side = max(EW,SN) 
mid_x, mid_y = LL[0]+EW/2.0, LL[1]+SN/2.0  #center location 
 
extent = [mid_x-side/2.0, mid_x+side/2.0, mid_y-side/2.0, mid_y+side/2.0]   # map coordinates, meters 
 
ax.set_extent(extent, crs=crg) 
ax.add_image(request, 8) 
 
# add a marker at center of the map 
plt.plot(mid_x, mid_y, marker='o', \ 
         color='red', markersize=10, \ 
         alpha=0.7, transform = crg) 
 
plt.show()