matlab：如何将ERA5(.nc)批量输出为Geotiff文件？

01 文件格式要求

文件样式形如2m_dewpoint_temperature_2000_01.nc

02 源码

1. clc
2. clear
3. %% 准备
4. in_dir = 'G:/ERA5/2m_temperature';  % 存储待处理nc文件的路径
5. out_dir='E:\MyTEMP\ERA5';  % 输出路径
6. start_year = 2000;  % 起始年份
7. end_year = 2000;  % 终止年份
8. var_wildcard_name = '2m_temperature';
9. var_name = 't2m';  % 变量名称
10. [status, msg] = mkdir(out_dir);  % 创建目录(存在则不创建)
11. srs=georefcells([-90 90],[-180 180],[1801 3600],'ColumnsStartFrom','north','RowsStartFrom','west');  % 坐标系信息
12. %{
13. georefcells函数用于针对一个规则栅格矩阵(确定的经纬度范围, 具体的行列数),创建一个
14. 地理参考对象.
15. 'ColumnsStartFrom','north' ==> 默认列方向是自北向南
16. 'RowsStartFrom','west'  ==> 默认行方向是自西向东
17. %}
18. %% 迭代获取变量的栅格矩阵并输出未Geotiff文件
19. % 检索nc文件
20. nc_paths = dir(fullfile(in_dir, '*.nc'));
21. for cur_year = start_year:end_year
22.     for cur_month = 1:1
23.         mdays = eomday(cur_year, cur_month);  % 获取当前年当前月份的天数
24.         % 读取当前循环nc文件的变量数组和基本信息
25.         cur_nc_path = fullfile(in_dir, sprintf('%s_%d_%02d.nc', var_wildcard_name, cur_year, cur_month));
26.         vinfo = ncinfo(cur_nc_path, var_name);
27.         [col, row, time] = deal(vinfo.Size(1), vinfo.Size(2), vinfo.Size(3));
28.         time_step = 0;  % 用于时间步计数
29.         for cur_day = 1:mdays
30.             for cur_hour = 1:24
31.                 time_step = time_step + 1;
32.                 var = ncread(cur_nc_path, var_name, [1, 1, time_step], [col, row, 1]);
33.                 %{
34.                 ncread(nc文件路径, 变量名称, 开始位置, 读取元素数)
35.                 [1, 1, time_step], 第一个维度(col)从1开始, 第二个维度(row)也是,
36.                 第三个维度(time)从第time_step开始
37.                 [col, row, 1], 第一个维度读取的元素数为col, 第二个为row, 第三个说明只读取一个时间点
38.                 %}
39.                 var = var';  % 转置-行列互换
40.                 var = [var(:, 1801:3600), var(:, 1:1800)];
41.                 %{
42.                 由于经度是0~360,这对应于地理坐标系就是0~180,-180~0
43.                 为了保证这与srs的地理参数设置一致(左上角点为-180°W,90°N)
44.                 %}
45.                 % 输出
46.                 cur_out_filename = sprintf('era5_%s_%d%02d%02d_%02d.tif', var_wildcard_name, cur_year, cur_month, cur_day, cur_hour);
47.                 cur_out_path = fullfile(out_dir, cur_out_filename);
48.                 geotiffwrite(cur_out_path, var, srs);
49.                 fprintf('当前处理: %s\n', cur_out_filename);
50.             end
51.         end
52.     end
53. end

复制代码

03 部分函数说明

3.1 georefcells函数

georefcells函数会返回一个针对地理坐标中规则像元栅格的默认参考对象。一般与geotiffwrite函数搭配利用(返回值可作为其函数的R参数传入)。
有三种传参方式：

• srs = georefcells(latlim,lonlim,rasterSize)
  

   latlim: 表现纬度范围,比方[-90, 90]表现从90°S~90°N;
    lonlim: 表现经度范围,比方[0, 180]表现从0°~180°E;
    rasterSize: 表现栅格大小,这里指代栅格图像的行列数,比方[1800, 3600]表现1800行3600列;
  

• srs = georefcells(latlim,lonlim,latcellextent,loncellextent)
  

这里的传参前面latlim,lonlim是一样,但是背面传入参数不一样。
   latcellextent: 表现纬度方向(y方向)上单个像元/栅格/单位的表现的范围，实际上就是纬度方向上的像元分辨率;
    loncellextent: 表现经度方向(x方向)上单个像元/栅格/单位的表现的范围，实际上就是纬度方向上的像元分辨率;譬如GPM IMERG降水影像的空间分辨率为0.1°×0.1°(经纬度方向上分辨率一致,范围全球),那么这里的就可以传入:
    srs = georefcells([-90, 90], [-180, 180],0.1,0.1)
  

• srs = georefcells(latlim,lonlim,___,Name,Value)
  

除了传入上述必要参数外，尚有一些可选参数通过关键字传参,这里包罗:

1. R = georefcells([-90, 90], [-180, 180], 0.1, 0.1,'ColumnsStartFrom','north')
2. R = georefcells([-90, 90], [-180, 180], 0.1, 0.1,'RowsStartFrom','east')

复制代码

  ColumnsStartFrom表现列方向上是从那里开始,默认是从北开始,即从北到南,如果栅格矩阵默认最上方是北边那么此处可以不举行关键字传参因为默认上方就是北(north);
    RowsStartFrom表现行方向上是从那里开始,默认是从西开始,即从西到东,如果栅格矩阵默认最左侧是东边那么此处可以不举行关键字传参因为默认左侧就是东(east);
  一般环境下，对于目前我所碰到的卫星遥感影像，大部分都是上北下南,左西右东，以是对于RowsStartFrom是必要指定为west.但是也有例外,对于MODIS的部分产品比方OMI/SWATH等产品部分存在南北颠倒即栅格矩阵上方是南极区域下方是北极区域.不外对于这种我一般还是不去更改这个ColumnsStartFrom参数,而是将读取的栅格矩阵举行南北极调换颠倒过来再举行后续操作.
3.2 eomday函数

很简单的一个函数，就是依据传入的年和月盘算得到当前年份当前月份下总共有多少天，官方说明即是：返回年份 Y 中指定月份 M 的末了一天(End of month day, eomday);
eomday传参如下：

1. E = eomday(Y,M)

复制代码

  Y: 表现年份;
    M: 表现月份;
  3.3 ncinfo函数

获取当前nc文件大概当前nc文件下指定变量/group的基本信息.
有如下三种传参方式：

1. finfo = ncinfo(source)
2. vinfo = ncinfo(source,varname)
3. ginfo = ncinfo(source,groupname)

复制代码

  source：表现nc文件的路径，不指定varname和group参数即返回该nc文件的基本信息;
    varname: 表现待查询的变量名称,注意此处变量的名称实在应当是变量在nc文件内的路径,比方t2m变量在group1/group2/t2m下，那么应该传入该路径而不是仅仅传入t2m这个变量名称;
    group_name: 表现待查询的group名称,说明与varname基本相似;
  3.4 ncread函数

获取当前nc文件指定变量名称的数组.
传参方式主要有下面两种：

1. vardata = ncread(source,varname)
2. vardata = ncread(source,varname,start,count)

复制代码

区别在于:如果你只是单纯想提取出整个变量数组，直接利用第一个传参方式即可;如果你必要提取变量数组中指定区域的部分数组出来，可以利用第二个传参方式(一般是在整个变量数组无法全部加载进电脑内存中，大概你必要针对变量数组中差别部分的小数组举行差别操作比方对于一个shape=(行,列,时间)的数组,大概对于差别时间的栅格矩阵必要单独处置处罚亦大概只必要某一个时间步的栅格矩阵)
   source: 表现nc文件的路径;
    varname: 表现待读取变量的名称,实际是变量在nc文件中的路径;
    start: 表现各个维度上的起始索引,比方对于shape=(行=100,列=200,波段数=3)的三维数组, 我们想要读取第二个波段的栅格矩阵可以设置[1, 1, 2]，再配合count参数(下面马上提及)设置为[100, 200, 1]即可;
    count: 表现沿每个维度读取的元素数;比方上述例子中必要提取第二个波段的栅格矩阵即从start中第一个维度(行)的索引1(matlab索引从1开始而非0)开始以后读取100个元素,类似地从列维度的索引1开始以后读取200个元素,从波段维度的索引2开始读取一个元素(即只读取一个元素就是第二个波段);
  3.5 geotiffwrite函数

传参方式: geotiffwrite(filename,A,R)
   filename: 输出Geotiff文件的输出路径;
    A: 输出的栅格矩阵,对于灰度图像大概单波段影像要求shape=(行,列);而对于RGB、多波段大概高光谱影像,要求shape=(行,列,波段);
    R: 地理参考对象,可以通过前面的georefcells函数构建得到;

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。