【Python基础绘图】Cartopy可视化Sentinel卫星影像

【Python基础绘图】Cartopy可视化Sentinel卫星影像

01 引言：

最近需要大量可视化哨兵影像，故研究一下，借助python cartopy批量化实现该需求，记录在此，分享给有需要的同学。
影像是通过GEE下载的【Harmonized Sentinel-2 MSI: MultiSpectral Instrument, Level-2A】的【B2、B3、B4】波段。

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 lib
import time
import cartopy.crs as ccrs
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
from osgeo import gdal


def cm2inch(x,y): 
    return x/2.54,y/2.54

t1 = time.time()

size1 = 10.5
fontdict = {'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.RasterXSize
rows = ds.RasterYSize
xmax = xmin+xres*cols
ymin = ymax+yres*rows
extent=[xmin,xmax,ymin,ymax]

data = ds.ReadAsArray().transpose((1, 2, 0))/3000

proj=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展示结果如下：

以上就是本期推文的全部内容了，如果对你有帮助的话，请‘点赞’、‘收藏’，‘关注’，你们的支持是我更新的动力。