Python + GDAL处理数据(2): 数据转存为栅格

接上篇
Python + GDAL 处理数据(1): NC文件的读取
老规矩，一上来先声明注释和导入包

# -*- coding:utf-8 -*-
from osgeo import gdal,osr,ogr,gdalconst
import os,datetime,random

# 读取数据文件
PrRCP85 = nc.Dataset("./pr_Amon_CanESM5_rcp85-cmip5_r5i1p1f1_gn_200601-210012.nc")

Year = 2015 
locals()["Pr"+str(Year)] = PrRCP85.variables['pr'][1,:,:] #这里是为了批量生成变量名，就不改了
locals()["Pr"+str(Year)] = np.asarray(locals()["Pr"+str(Year)])
locals()["Pr"+str(Year)][np.where(locals()["Pr"+str(Year)] == 1e+20)] = np.nan #异常值处理
locals()["Pr"+str(Year)] = locals()["Pr"+str(Year)] * 1000 * 100 * 100
Lat = PrRCP85.variables['lat'][:]
Lon = PrRCP85.variables['lon'][:]

PrOutTif = "./Tiff Store/RCP85/Raw Data Tiff Store/ACCESS-ESM1-5/MAP/Pr_%s.tif"%Year
#Tif影像的信息
LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()] 
xsize = (LonMax-LonMin)/len(Lon)
ysize = (LatMax-LatMin)/len(Lat)

#构建栅格
locals()["Pr"+str(Year)+"_ds"] = gdal.GetDriverByName('Gtiff').Create(PrOutTif,len(Lon),len(Lat),1,gdal.GDT_Float32) 
locals()["Pr"+str(Year)+"_ds"].SetGeoTransform((LonMin,xsize, 0, LatMin, 0, ysize))

#设置地理坐标系和投影坐标系，这里设置为WGS_84,索引号4326
srs = osr.SpatialReference()# 获取地理坐标系统信息，用于选取需要的地理坐标系统
srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
locals()["Pr"+str(YearsNum[Ind+80])+"_ds"].SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息
locals()["Pr"+str(YearsNum[Ind+80])+"_ds"].GetRasterBand(1).SetNoDataValue(-9999) #设置空值玮-9999
locals()["Pr"+str(YearsNum[Ind+80])+"_ds"].GetRasterBand(1).WriteArray(locals()["Pr"+str(YearsNum[Ind+80])]) # 将数据写入内存
locals()["Pr"+str(YearsNum[Ind+80])+"_ds"].GetRasterBand(1).GetStatistics(0,1)# 添加统计信息
locals()["Pr"+str(YearsNum[Ind+80])+"_ds"].FlushCache() #将数据写入硬盘
del locals()["Pr"+str(YearsNum[Ind+80])+"_ds"] #关闭栅格文件，解除文件占用

到这里，就完成了数据绘制栅格部分