全球干旱指数PDSI可从以下网站下载：
http://climexp.knmi.nl/selectfield_obs2.cgi?id=someone@somewhere
image.png
从该网站上可直接下载全球1901-2016年的逐月PDSI数据，文件为nc格式，为进一步和其他栅格数据进行计算，需要将nc文件转变为tif文件，因此本文提供一种能够批量转换的处理方式。
首先准备个样例数据
ncdisp('H:\Global\PDSI\scPDSI.cru.3.25.bams2017.GLOBAL.1901.2016.nc');
data1=ncread('H:\Global\PDSI\scPDSI.cru.3.25.bams2017.GLOBAL.1901.2016.nc','scpdsi');
data3=data1(:,:,1);
data4=rot90(data3);
data5=flipud(data4);
data5(isnan(data5))=-999;
dlmwrite('样例1.txt',data5,'\t',1,1);
经过上述转换后可得到文本格式的样例数据，结果如下：
image.png
然后在文本中他添加经纬度信息，结果如下所示：
image.png
添加经纬度信息后，利用arcgis的ASCII码转raster功能，将该文本转变为栅格文件，并进一步输出为tif，假设文件名为样例.tif。
然后加载带有投影信息的其他栅格文件，对样例.tif文件定义投影，投影方式与加载进来的栅格文件保存一致。
经过上述步骤可得到样例数据。接下来进行批量的转换。
[aaaaa,R]=geotiffread('H:\Global\PDSI\example1.tif');%先导入纬度数据
info=geotiffinfo('H:\Global\PDSI\example1.tif');
data=ncread('H:\Global\PDSI\scPDSI.cru.3.25.bams2017.GLOBAL.1901.2016.nc','scpdsi');
for year=1901:2016
data1=data(:,:,1+12*(year-1901):12*(year-1900)); %得到每年的12个月数据
data3=sum(data1,3)/12;
data4=rot90(data3);
data5=flipud(data4);
filename=strcat('H:\Global\PDSI\年尺度pdsi\全球',int2str(year),'年PDSI.tif');
geotiffwrite(filename,data5,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
for mon=1:12
data2=data1(:,:,mon);
data4=rot90(data2);
data5=flipud(data4);
filename=strcat('H:\Global\PDSI\月尺度的pdsi\全球',int2str(year),'_',int2str(mon),'月PDSI.tif');
geotiffwrite(filename,data5,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
end
end
通过上述代码即可实现批量处理nc格式的PDSI文件了, 在新的版本中，可以不用样例文件直接进行提取，在VX公众号地学分析与算法中，结果如下：
image.png
image.png
CRU文件的处理过程
data=ncread('H:\Global\CRU\cru_ts4.01.1901.2016.tmx.dat.nc','tmx');
data3=data(:,:,1);
data4=rot90(data3);
data5(isnan(data5))=-999;
dlmwrite('样例2.txt',data5,'\t',1,1);
其余步骤同上构建一个有投影的样例tif数据，批量读取与转换
[aaaaa,R]=geotiffread('H:\Global\CRU\样例2.tif');%先导入纬度数据
info=geotiffinfo('H:\Global\CRU\样例2.tif');
data=ncread('H:\Global\CRU\cru_ts4.01.1901.2016.tmx.dat.nc','tmx');
for year=1901:2016
data1=data(:,:,1+12*(year-1901):12*(year-1900)); %得到每年的12个月数据
data3=sum(data1,3)/12; %对年数据求平均值，得到年平均最大气温，如果是降水，则直接去掉/12
data4=rot90(data3);
filename=strcat('H:\Global\CRU\tif\year\CRU',int2str(year),'_Tmx.tif');
geotiffwrite(filename,data4,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
for mon=1:12
data2=data1(:,:,mon);
data4=rot90(data2);
filename=strcat('H:\Global\CRU\tif\month\CRU',int2str(year),'_',int2str(mon),'_Tmax.tif');
geotiffwrite(filename,data4,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
end
end
更多需求，请查看个人介绍
