01 引言：
02 代码如下：
# -*- encoding: utf-8 -*-'''@File : demo2.py@Time : 2023/02/12 13:57:10@Author : HMX@Version : 1.0@Contact : kzdhb8023@163.com'''# here put the import libimport timeimport cartopy.crs as ccrsimport matplotlib as mplimport matplotlib.pyplot as pltimport numpy as npfrom cartopy.io.shapereader import Readerfrom cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatterfrom osgeo import gdaldef cm2inch(x,y):return x/2.54,y/2.54t1 = time.time()size1 = 10.5fontdict = {'weight': 'bold','size':size1,'family':'SimHei'}mpl.rcParams.update({'text.usetex': False,'font.family': 'stixgeneral','mathtext.fontset': 'stix',"font.family":'serif',"font.size": size1,"mathtext.fontset":'stix',"font.serif": ['Times New Roman'],})shp_path = r'E:\CODE\study\image_plot\data\奎文.shp'tif_path = r'E:\CODE\study\image_plot\data\奎文区2022-08.tif'ds = gdal.Open(tif_path)geotransform = ds.GetGeoTransform()xmin = geotransform[0]ymax = geotransform[3]xres = geotransform[1]yres = geotransform[-1]cols = ds.RasterXSizerows = ds.RasterYSizexmax = xmin+xres*colsymin = ymax+yres*rowsextent=[xmin,xmax,ymin,ymax]data = ds.ReadAsArray().transpose((1, 2, 0))/3000proj=ccrs.PlateCarree()fig,ax1 = plt.subplots(1, 1,figsize = (8,6),dpi=100, subplot_kw={'projection': proj})# 绘制影像ax1.imshow(data.clip(0,1), extent=extent)# 绘制矢量边界ax1.add_geometries(Reader(shp_path).geometries(), proj, edgecolor='r', facecolor='none',alpha=1, linewidth=0.65)ax1.set_xticks([np.round(i,2) for i in np.arange(extent[0], extent[1]+xres, 0.05)], crs = proj)ax1.set_yticks([np.round(i,2) for i in np.arange(extent[2], extent[3]+yres,0.05)], crs = proj)ax1.xaxis.set_major_formatter(LongitudeFormatter())ax1.yaxis.set_major_formatter(LatitudeFormatter())ax1.minorticks_on()plt.tight_layout()plt.savefig(r'E:\CODE\study\image_plot\data\奎文.png',dpi = 600)t2 = time.time()print('共计用时{:.2f}s'.format(t2-t1))plt.show()
03 结果如下：
对比arcgis展示结果如下：
戳这里关注我
请点赞、在看、关注，你们的支持是我更新的动力。
