同学写的代码,学习一下
import numpy as np
from osgeo import gdal,osr
import rasterio.plot
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import matplotlib as mpl
import matplotlib.ticker as mticker
from matplotlib_scalebar.scalebar import ScaleBar
from rasterio.plot import show
import north
from geopy import distance
import os
from north import add_north
import seaborn as sns
def getSRSPair(dataset):
'''
获得给定数据的投影参考系和地理参考系
:param dataset: GDAL地理数据
:return: 投影参考系和地理参考系
'''
prosrs = osr.SpatialReference()
prosrs.ImportFromWkt(dataset.GetProjection())
geosrs = prosrs.CloneGeogCS()
return prosrs, geosrs
def geo2lonlat(dataset, x, y):
'''
将投影坐标转为经纬度坐标(具体的投影坐标系由给定数据确定)
:param dataset: GDAL地理数据
:param x: 投影坐标x
:param y: 投影坐标y
:return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)
PROJCS["WGS 84 / UTM zone 45N"]
'''
prosrs, geosrs = getSRSPair(dataset)
ct = osr.CoordinateTransformation(prosrs, geosrs)
coords = ct.TransformPoint(x, y)
return coords[:2]
# 创建画布
fig = plt.figure(figsize=(15,10))
# fig = plt.figure(figsize=(6,4))
image_names = os.listdir('D:\\影像练习\\NDVI输出')
#坐标轴范围
ds = gdal.Open('D:\\影像练习\\001.tif')
gt = ds.GetGeoTransform()
print(gt)
coordss = geo2lonlat(ds, gt[0], gt[3])
coordsss = geo2lonlat(ds, gt[0] + ds.RasterXSize * gt[1], gt[3] + ds.RasterYSize * gt[5])
extent = (coordss[1], coordsss[1],coordsss[0], coordss[0])
# 生成三个子画布,并生成标题、比例尺、网格
dx = distance.great_circle((extent[0],extent[2]),(extent[0]+1,extent[2])).km
axs = []
for i in range(1,4):
axs.append(fig.add_subplot(int('13'+str(i))))
axs[i-1].set_title(image_names[i-1][6:10]+'-'+image_names[i-1][10:12]+'-'+image_names[i-1][12:14]+'\n', size=15)
scalebar = ScaleBar(dx=dx,units='km',box_alpha=0,scale_loc='bottom',location='lower right',color='k')
plt.gca().add_artist(scalebar)
plt.grid(linestyle='--',color='gray',alpha=0.5)
# 定义NDVI极值
vmin = 0.001965923984272608
vmax = 0.35251450676982593
# vmin = 0.0
# vmax = 1.0
# 颜色空间和值空间的映射
norm = mpl.colors.Normalize(vmin=vmin,vmax=vmax)
cmap = ListedColormap(['#FFF8DC','#008B8B','PeachPuff','Tomato','blue'])
# cmap = 'gray'
# 读取之前生成的NDVI的tif数据
ims = []
rasters = []
for image_name in image_names:
im = rasterio.open('D:\\影像练习\\NDVI输出'+'\\'+image_name)
ims.append(im)
raster = im.read()
for i in range(len(raster[0])):
for j in range(len(raster[0][i])):
# print(raster[0][i][j])
if vmax-(vmax-vmin)/4<=raster[0][i][j]<vmax:
raster[0][i][j] = ((vmax-(vmax-vmin)/4)+(vmax))/2
elif vmax-(vmax-vmin)/2<=raster[0][i][j]<vmax-(vmax-vmin)/4:
raster[0][i][j] = ((vmax-(vmax-vmin)/2)+(vmax-(vmax-vmin)/4))/2
elif vmax-3*(vmax-vmin)/4<=raster[0][i][j]<vmax-(vmax-vmin)/2:
raster[0][i][j] = ((vmax-(vmax-vmin)/2)+(vmax-3*(vmax-vmin)/4))/2
elif vmin<=raster[0][i][j]<vmax-3*(vmax-vmin)/4:
raster[0][i][j] = ((vmin)+(vmax-3*(vmax-vmin)/4))/2
# print(raster)
rasters.append(raster)
# 在各子图上添加NDVI地图并绘制指北针
raster_num = 0
transform = (0.00016, 0.0, extent[0],0.0, -0.00016, extent[3])
print(extent)
print(ims[0].profile)
for ax in axs:
show(rasters[raster_num],ax=ax,cmap=cmap,norm=norm,extent=extent,transform=transform)
ax.yaxis.set_major_formatter(mticker.FormatStrFormatter('%0.2f°N'))
ax.xaxis.set_major_formatter(mticker.FormatStrFormatter('%0.2f°E'))
raster_num += 1
if raster_num > 2:
raster_num = 0
add_north(ax) #添加指北针
#添加图例
pos_size = fig.add_axes([0.17,0.15,0.7,0.03]) #图例的位置和尺寸
cb1 = plt.colorbar(mpl.cm.ScalarMappable(norm=norm,cmap=cmap),orientation='horizontal',label='NDVI',cax=pos_size)
# cb1.set_ticks(np.arange(vmax,vmin,-0.25)) #设置图例刻度
cb1.set_ticks(np.arange(vmax,vmin,-0.07)) #设置图例刻度
#显示并保存专题地图
plt.savefig('D:/影像练习/渲染/NDVI_lon_lat.jpg')
plt.savefig('D:/影像练习/渲染/NDVI_lon_lat.pdf')
plt.show()