python批量读取中国降水格点数据转tif(附数据)
中国地面降水月值0.5°×0.5°格点数据集(V2.0)196101-201305
这个数据提供了1961-2013年的长时期站点观测,由气象站插值得到,质量较高。
官方的描述是:
中国地面降水月值0.5°×0.5°格点数据集文件命名由数据集代码(SURF_CLI_CHN_PRE_MON_GRID_0.5)、年份、月份标识(YYYYMM)组成。 时间范围:1961年1月-最新
时间分辨率:逐月 3.4.空间属性 3.4.1.地理范围 地理范围描述:中国
最西经度:72°E
最东经度:136°E
最北纬度:54°N
最南纬度:18°N
这个数据我也免费提供给大家,见文末。
然而数据是这样的:
我们在excel尝试打开:
天天被国产数据搞破防,关键是这个还不是固定空格分隔的数据,而且txt还多出没用的几行,给读取带来麻烦:
直接上代码:
import pandas as pd
import numpy as np
import os
dir = 'D:/Acdemic/xibei/data_1/grid_prcp/'
txtLists = os.listdir(dir)
files = list(filter(lambda x: x[-4:] in ['.txt'], txtLists))
df = pd.DataFrame()
i = 0;
prcp = np.zeros((72, 128, len(files)))
for file in files:
data = pd.read_table(dir+file, sep='\s+', header=None, index_col=False, skiprows=6)
print(str(i) + ': ' + file)
prcp[:, :, i] = np.array(data)
i += 1
这样会读取文件下的所有数据,转为ndarray,并print出序号和文件名,结果为:
0: SURF_CLI_CHN_PRE_MON_GRID_0.5-200001.txt
1: SURF_CLI_CHN_PRE_MON_GRID_0.5-200002.txt
2: SURF_CLI_CHN_PRE_MON_GRID_0.5-200003.txt
3: SURF_CLI_CHN_PRE_MON_GRID_0.5-200004.txt
4: SURF_CLI_CHN_PRE_MON_GRID_0.5-200005.txt
5: SURF_CLI_CHN_PRE_MON_GRID_0.5-200006.txt
6: SURF_CLI_CHN_PRE_MON_GRID_0.5-200007.txt
7: SURF_CLI_CHN_PRE_MON_GRID_0.5-200008.txt
8: SURF_CLI_CHN_PRE_MON_GRID_0.5-200009.txt
9: SURF_CLI_CHN_PRE_MON_GRID_0.5-200010.txt
10: SURF_CLI_CHN_PRE_MON_GRID_0.5-200011.txt
11: SURF_CLI_CHN_PRE_MON_GRID_0.5-200012.txt
接下来求所有降水数据平均值,由于我想求年总降水,数据是月均值。乘12即可。
prcp1 = np.mean(prcp, 2) * 12
prcp1[prcp1 < -1000] = None
最后将其导为TIF,根据文档的属性信息,手动记录下左下角和右上角坐标和分辨率:
LonMin = 72.5
Lon_Res = 0.5
Lat_Res = 0.5
LatMax = 54.5
driver = gdal.GetDriverByName('GTiff')
out_tif_name = 'D:/Acdemic/xibei/data_1/grid_mean.tif'
out_tif = driver.Create(out_tif_name,N_Lon,N_Lat,1,gdal.GDT_Float32)# 设置影像的显示范围#-Lat_Res一定要是-的
geotransform = (LonMin,Lon_Res, 0, LatMax, 0, -Lat_Res)
out_tif.SetGeoTransform(geotransform)#获取地理坐标系统信息,用于选取需要的地理坐标系统
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)# 定义输出的坐标系为"WGS 84",AUTHORITY["EPSG","4326"]
out_tif.SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息#数据写出
out_tif.GetRasterBand(1).WriteArray(prcp1) # 将数据写入内存,此时没有写入硬盘
out_tif.FlushCache() # 将数据写入硬盘
out_tif = None # 注意必须关闭tif文件
可以看到已经导出成功了:
现在我把这个数据共享:
现在我把这个数据共享:
后台回复【中国降水格点数据集】