python 渲染

   同学写的代码，学习一下

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()