ax.add_geometries无法加载shp文件

作者在撰写毕业论文时遇到cartopy包中shp文件在地图上不显示的问题,经排查发现是由于绘图范围与shp文件的投影范围不匹配。通过使用Arcgis的project工具将shp文件投影转换至地理坐标系,解决了这个问题。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

最近在写毕业论文,上手练习一下cartopy包。

问题:

正常运行ax.add_geometries代码,但是shp文件在结果图中一直不显示。一开始还以为是投影有问题,代码中设置的是简易圆柱(Plate Carrée) 投影,但后面也尝试了其他投影,依然显示不了。

ax.add_geometries(Reader(shp_path).geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=0.65)

运行结果如下:

问题发现:

原来是子图中的绘图范围与Shp的数值范围不匹配!

即:我之前提供的shp文件是投影过的,检查代码如下:

import geopandas as gpd
# 读取Shapefile
gdf = gpd.read_file(shp_path)
# 获取边界
bounds = gdf.total_bounds
print("Shapefile bounds:", bounds)
#结果: [11448677.0 3047834.8 11625463.8 3218060.9 ]

然而,我的绘图范围确是用经纬度表示的,检查代码如下:

xlim = ax.get_xlim()#这个ax是我子图的轴
ylim = ax.get_ylim()
print("Current map view limits:", xlim, ylim)
#结果:Current map view limits: (102.5, 105.0) (27.0, 29.5)

解决:

将shp文件投影更替为地理坐标系即可,我是用的Arcgis的project工具。然后,再检查shp文件的数值范围时,其范围就落入绘图范围中啦!

#同样的代码
import geopandas as gpd
# 读取Shapefile
gdf = gpd.read_file(shp_path)
# 获取边界
bounds = gdf.total_bounds
print("Shapefile bounds:", bounds)
#结果更新: [102.9 27.5 104.7 28.9]

绘图结果:

完结撒花!

对代码修改1、色带采用地形渐变色带(Terrain)terrain(Matplotlib内置):低海拔用绿色 → 中海拔黄棕色 → 高海拔灰白色2、对主图和子图的大小进行适当调整,主图可以小一些,子图可以大一些,但是要中间对齐3、色带匹配在主图的下方并且与主图的距离控制位0.02,且长度匹配。4、要求最后输出的图片主图与子图的图幅一定要完整。具体代码为:import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.io.shapereader import Reader from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter import xarray as xr import geopandas as gpd import matplotlib.path as mpath from matplotlib import rcParams # 配置字体和大小 config = {"font.family": 'Times New Roman', "font.size": 16, "mathtext.fontset": 'stix'} rcParams.update(config) # 定义区域范围 region = [70, 140, 10, 60] # 打开数据集 with xr.open_dataset(r'C:\Users\asus\Desktop\ETOPO2v2c_f4.nc') as ds: ds = ds.sel(x=slice(60, 150), y=slice(0, 60)) lon, lat = np.meshgrid(ds['x'], ds['y']) dem = ds['z'].to_numpy() # 创建画布和投影 fig = plt.figure(dpi=600, figsize=(16, 9)) proj = ccrs.PlateCarree() # 第一个子图(主图) ax = fig.add_axes([0.05, 0.05, 0.4, 0.85], projection=proj) # 调整子图位置和大小 ax.set_extent(region, crs=proj) ax.set_xticks(np.arange(region[0], region[1] + 1, 10), crs=proj) # 主图经度刻度 ax.set_yticks(np.arange(region[2], region[3] + 1, 10), crs=proj) # 主图纬度刻度 ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False)) ax.yaxis.set_major_formatter(LatitudeFormatter()) # 添加地图要素 ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, zorder=2, color='k') ax.add_feature(cfeature.LAKES.with_scale('50m')) ax.add_feature(cfeature.OCEAN.with_scale('50m')) ax.add_feature(cfeature.LAND.with_scale('50m')) # 添加九段线 ax.add_geometries(Reader(r'C:\Users\asus\Desktop\中国省级地图GS(2019)1719号\九段线GS(2019)1719号.shp').geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor='k', linewidth=1) # 绘制 DEM 数据,使用更符合高程的色调 cmap = plt.cm.terrain # 使用 terrain 色调 c11 = ax.pcolormesh(lon, lat, dem, transform=ccrs.PlateCarree(), cmap=cmap, vmin=0, vmax=6000) # 读取中国省级地图并裁剪 gdf = gpd.read_file(r'C:\Users\asu
03-12
f_path = r"E:\gra_thesis\sum_pre_data_new\grid_nc\AMJ_pre_total_precip.nc" f = xr.open_dataset(f_path) f # %% lon = f['lon'] lat = f['lat'] data= f['precip'] data_mean = np.mean(data, 0) # %% shp_path = r"C:\Users\86133\Desktop\thesis\2020国家级行政边界\China_province.shp" sf = shapefile.Reader(shp_path) shp_reader = Reader(shp_path) sf.records() region_list = [110000, 120000, 130000,140000,150000,210000,220000, 230000, 310000, 320000,330000,340000,350000,360000, 370000, 410000, 420000,430000,440000,450000,460000, 500000, 510000, 520000,530000,540000,610000,620000, 630000, 640000, 650000,710000,810000,820000] # %% proj = ccrs.PlateCarree() extent = [105, 125, 15, 30] fig, ax = plt.subplots(1, 1, subplot_kw={'projection': proj}) ax.set_extent(extent, proj) # ax.add_feature(cfeature.LAND, fc='0.8', zorder=1) ax.add_feature(cfeature.COASTLINE, lw=1, ec="k", zorder=2) ax.add_feature(cfeature.OCEAN, fc='white', zorder=2) ax.add_geometries(shp_reader.geometries(), fc="None", ec="k", lw=1, crs=proj, zorder=2) ax.spines['geo'].set_linewidth(0.8) ax.tick_params(axis='both',which='major',labelsize=9, direction='out',length=2.5,width=0.8,pad=1.5, bottom=True, left=True) ax.tick_params(axis='both',which='minor',direction='out',width=0.5,bottom=True,left=True) ax.set_xticks(np.arange(105, 130, 5)) ax.set_yticks(np.arange(15, 40, 5)) ax.xaxis.set_major_formatter(LongitudeFormatter()) ax.yaxis.set_major_formatter(LatitudeFormatter()) cf = ax.contourf(lon, lat, data_mean, extend='both', cmap='RdBu') cb = fig.colorbar(cf, shrink=0.9, pad=0.05)解释这段代码
04-20
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值