ax.add_geometries无法加载shp文件

最近在写毕业论文,上手练习一下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]

绘图结果:

完结撒花!

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、付费专栏及课程。

余额充值