全球干旱指数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文件了。结果如下：

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