matlab指数转换,基于Matlab的干旱指数PDSI及CRU全球气象nc数据的处理过程-由nc文件转变为tif文件...

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

更多需求，请查看个人介绍