CRU气候数据处理

记录一下自己的数据处理流程，以供此后翻阅。

一、CRU数据下载

网址：High-resolution gridded datasets (uea.ac.uk)https://crudata.uea.ac.uk/cru/data/hrg/下载流程参考：1901年-2020年全球气象数据CRU TS 介绍、下载与使用教程 - 知乎 (zhihu.com)

下载文件（以pre为例）：cru_ts4.07.1901.2022.pre.dat.nc.gz

二、CRU数据处理

多年nc文件批量转为逐月tif文件 

代码来自：matlab批量读取nc文件并转为tif_哔哩哔哩_bilibili

clc;
clear;

%读取nc文件基本信息
inpath = 'E:\variables\CRU\PRE\';  %nc文件所在文件夹
info1 = ncinfo('E:\variables\CRU\PRE\cru_ts4.07.1901.2022.pre.dat.nc');  %nc文件基本信息

%先在arcgis中读取nc文件中的一个时间点数据，并导出为tif，以便后面为nc文件批量导出提供坐标系
maskname = 'D:\Arcgis\tif\cru_pre_0.5×0.5.tif';
[A,R] = geotiffread(maskname);  %主要目的是获取R，也就是坐标系
info = geotiffinfo(maskname);

%获取nc文件基本信息
infile = strcat(inpath,'cru_ts4.07.1901.2022.pre.dat.nc');
lon = ncread(infile,'lon');
lat = ncread(infile,'lat');
D = ncread(infile,'pre');   %读取变量D的内容

%按时间分别读取
for year = 2000:2020  %年循环 
    for month = 1:12  %月循环
        n = (year-1901)*12+month; %1901为数据起始年份
        data = D(:,:,n);  %逐月读取
        data = rot90(data);  %rot90逆时针旋转90°/flipud上下翻转/fliplr左右翻转  
                             %有时候数据行列数反了，需要旋转数据来调整
        export2tif = ['E:\variables\CRU\PRE\tif\pre_',int2str(year),'_',int2str(month),'_0.5×0.5.tif'];  %导出路径和名字
        geotiffwrite(export2tif,data,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);  %批量导出为tif
    end
    disp([year, month]);
end

中间出现过的问题：先在arcgis中导出的tif空间分辨率不能更改，否则报错。

给矩阵添加地理信息_r.rastersize is inconsistent with rastersize.-CSDN博客
 

补充：nc文件批量转为逐年tif文件

参考：使用matlab将多年温度数据在一个nc文件当中，按月和年分别转为tif文件_一个tif可以存储多个月份的数据吗-CSDN博客

pre  使用sum函数

clc;
clear;

% 读取 nc 文件基本信息
inpath = 'E:\variables\CRU\PRE\';  % nc 文件所在文件夹
ncFile = 'cru_ts4.07.1901.2022.pre.dat.nc';  % nc 文件路径
info = ncinfo(fullfile(inpath, ncFile));  % nc 文件基本信息

% 先在 ArcGIS 中读取 nc 文件中的一个时间点数据，并导出为 tif，以便后面为 nc 文件批量导出提供坐标系
maskname = 'D:\Arcgis\tif\cru_pre_0.5×0.5.tif';
[A, R] = geotiffread(maskname);  % 主要目的是获取 R，也就是坐标系
info = geotiffinfo(maskname);

% 获取 nc 文件基本信息
lon = ncread(fullfile(inpath, ncFile), 'lon');
lat = ncread(fullfile(inpath, ncFile), 'lat');
D = ncread(fullfile(inpath, ncFile), 'pre');   % 读取变量 D 的内容

% 按年分别读取
for year = 2000:2020  % 年循环
    startMonth = (year - 1901) * 12 + 1;
    endMonth = startMonth + 11;  % 对应年份的最后一个月份
    yearlyData = D(:, :, startMonth:endMonth);  % 提取整个年份的数据
    yearlyData = rot90(yearlyData);  %逆时针旋转图像，这里转出tif后可以先看一下数据的朝向以及是否需要翻转
    				                 %每种数据需要翻转和旋转的角度可能不一样，需要自己确定
    yearlyData(yearlyData == -32768) = NaN;  % 无数据区域设为 NaN

    % 计算年降水量累计值
    yearlySumData = sum(yearlyData, 3);
    
    % 保存每年累计值为.tif文件
    export2tif = ['E:\variables\CRU\PRE\yearly_tif\pre_', int2str(year), '_0.5×0.5.tif'];  % 导出路径和名字
    geotiffwrite(export2tif,yearlySumData,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);  % 批量导出为tif
    disp(year);
end

 tmp  使用mean函数  其余同pre

clc;
clear;

% 读取 nc 文件基本信息
inpath = 'E:\variables\CRU\TMP\';  % nc 文件所在文件夹
ncFile = 'cru_ts4.07.1901.2022.tmp.dat.nc';  % nc 文件路径
info = ncinfo(fullfile(inpath, ncFile));  % nc 文件基本信息

% 先在 ArcGIS 中读取 nc 文件中的一个时间点数据，并导出为 tif，以便后面为 nc 文件批量导出提供坐标系
maskname = 'D:\Arcgis\tif\cru_tmp_0.5×0.5.tif';
[A, R] = geotiffread(maskname);  % 主要目的是获取 R，也就是坐标系
info = geotiffinfo(maskname);

% 获取 nc 文件基本信息
lon = ncread(fullfile(inpath, ncFile), 'lon');
lat = ncread(fullfile(inpath, ncFile), 'lat');
D = ncread(fullfile(inpath, ncFile), 'tmp');   % 读取变量 D 的内容

% 按年分别读取
for year = 2000:2020  % 年循环
    startMonth = (year - 1901) * 12 + 1;
    endMonth = startMonth + 11;  % 对应年份的最后一个月份
    yearlyData = D(:, :, startMonth:endMonth);  % 提取整个年份的数据
    yearlyData = rot90(yearlyData);  %逆时针旋转图像，这里转出tif后可以先看一下数据的朝向以及是否需要翻转
    				                 %每种数据需要翻转和旋转的角度可能不一样，需要自己确定
    yearlyData(yearlyData == -32768) = NaN;  % 无数据区域设为 NaN

    % 计算年均温
    yearlyMeanData = mean(yearlyData, 3);
    
    % 保存每年平均值为.tif文件
    export2tif = ['E:\variables\CRU\TMP\yearly_tif\tmp_', int2str(year), '_0.5×0.5.tif'];  % 导出路径和名字
    geotiffwrite(export2tif,yearlyMeanData,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);  % 批量导出为tif
    disp(year);
end

二、CRU数据重采样

将得到的tif文件导入arcgis

利用模型构建器进行批量重采样至0.05°

参考：【ArcGIS教程】（113）模型构建器（9）——批量影像重采样 - 知乎 (zhihu.com)

利用arcgis的modelbuilder来批量对影像进行重采样_arcgis批量重分类-CSDN博客

三、提取站点CRU值

导入站点经纬度坐标，建立shp文件，而后多值提取至点。

遇到问题：shapefile 格式将字段名称的最大长度限制为 10 个字符。 因此，对于追加到输入 shapefile 属性表中的任何字段，其名称都将被截断并获得唯一值。 如果名称很长或很相似，则可能导致各字段间难以区分。 在这种情况下，建议您将输入 shapefile 复制到文件地理数据库，然后将要素类用于分析。

ArcGIS Help 10.2 - 多值提取至点 (空间分析)

解决：将shp文件导入地理数据库

https://jingyan.baidu.com/article/3065b3b6e883daffcef8a451.html

在站点的属性表即可看到pre信息

导出属性表：表转excel工具有限制（ 小于65535 行和 256 列），直接全选复制到excel表。