1 基于MATLAB的地理数据的读取

红木棉

中山大学自然地理学本硕

本章为专栏《基于MATLAB的地理数据处理》的内容，关于该专栏的目录初步安排请 戳此跳转 。

本章内容所涉及到的数据、数据说明与代码文件，请到百度网盘下载。本章列出的代码均在"A_load_data.m"里。

链接： https:// pan.baidu.com/s/198G2CK ull8Nwq052FrIR5g 提取码：mian

发布时间：2022年4月24日

修改时间：2022年5月3日(1.4节对lat_nc进行一次flipud，否则纬度朝向相反)

作者：棉

1.1 表格数据(excel，text)

在地理数据采样过程中，人工或自动站点采样的数据通常记录在excel、text以及csv等格式的数据中。MATLAB读取这类数据的函数为 readtable ( MATLAB官方文档 )，该函数读取上述格式的数据后，在工作区里所显示的数据类型都为table，因此本文统一称这类数据为表格数据。

首先，根据本文提供的代码，将MATLAB_Course文件夹放至自己计算机对应的位置，并将代码文件中的Path变量更改为自己所存放的路径。下面以存放在计算机的D盘为例：

clear;clc
% 指定MATLAB_Course的路径
Path = 'D:\MATLAB_Course';
cd(Path)

在 '\MATLAB_Course\Data\table_CMA_DATA\' 里，本文提供了2022年4月10日–4月15日珠三角部分国家级气象站点的气象观测数据，DATA_CMA_CHN_MUL_HOR.txt。这些气象站点的经纬度信息也分别保存在Location_Station.txt和Location_Station.xlsx里。

接下来，我们使用 readtable 函数读取DATA_CMA_CHN_MUL_HOR.txt：

data_txt = readtable('.\Data\table_CMA_DATA\DATA_CMA_CHN_MUL_HOR.txt');

读取后结果如图1所示，data_txt为3024×8的table类型数据。table类型数据有个表头，分别显示每列数据所代表的变量名称，依次为Station_ID、Year、Month、Day、Hour、PRS_hPa (气压/百帕)、TEM_deg (气温/摄氏度)、PRE_1h_mm (降水/毫米)。

接下来再尝试分别读取Location_Station.txt和Location_Station.xlsx。

loc_txt = readtable('.\Data\table_CMA_DATA\Location_Station.txt');
loc_xls = readtable('.\Data\table_CMA_DATA\Location_Station.xlsx');

关于如何进一步处理MATLAB的table类型数据，请大家阅读官方文档，本文仅提供一些浅层的提取处理。以读取后的loc_txt变量为例：

test1 = loc_txt{:,2:4};
test2 = loc_txt(:,2:4);
test3 = loc_txt.Station_Name;
test4 = loc_txt.Station_ID;

运行上述命令后的结果如图2所示：

具体请自行操作后查看test1、test2、test3和test4有何区别，一般来说，我们都是要提取出double类型的向量或矩阵才能方便地进行后续的处理与计算。大家还可以试试下面的会不会出错，思考下为什么？

test5 = loc_txt{:,1:2};
test6 = loc_txt(:,1:2);

1.2 矢量数据(shapefile)

经常使用Arcgis的人对shapefile这种文件格式应该并不陌生，这也就是我们常见的后缀为.shp的文件。在 '\MATLAB_Course\Data\shp_GuangZhou_Part\' 里，本文提供了广州市行政区的矢量边界数据(为了方便教学，我去掉了两个小岛屿的边界，所以严格意义上说这不是完整的广州市行政区边界)。

值得注意的一点是，在地理数据处理过程中，我们通常会把数据统一处理为地理坐标系，也就是经纬坐标系，而大多数矢量边界数据所采用的坐标系为投影坐标系，所以建议在用MATLAB读取之前先用Arcgis等软件将矢量数据转为WGS84经纬坐标系。坐标系转换后，即可用MATLAB的 shaperead 函数直接读取shp格式的数据：

gz_shp = shaperead('.\Data\shp_GuangZhou_Part\GuangZhou_Part.shp');

如图3所示，读取后的gz_shp变量的数据类型为struct类型，里面包含的X和Y即我们需要的坐标。

gz_shp.X % 经度
gz_shp.Y % 纬度

用上述两个命令即可提取出gz_shp矢量数据的坐标向量。

1.3 分层数据格式(HDF)

分层数据格式英文全称为Hierarchical Data Format，简称为HDF。该文件格式由美国国家超级计算应用中心(NCSA)开发，旨在帮助科学家更便捷地共享数据，而无需考虑底层的计算环境。美国国家宇航局(NASA)在此文件格式基础上进一步开发出HDF-EOS子集(EOS: Earth Observation System)，用于记录地理数据。NASA大多数地理数据，包括MODIS、SMAP、VIIRS、AMSR、OCO-2等均是采用HDF格式。因此，基于MATLAB的读取HDF数据的方法十分重要。本文分别以MODIS产品(HDF格式)与SMAP产品(HDF 5.0，简称h5)的读取为例进行介绍。

(1) 基于MATLAB的MODIS数据的读取(HDF)

在 '\MATLAB_Course\Data\hdf_MOD13C2_NDVI\' 里，本文提供了2021年7月的NDVI数据，所用产品为Vegetation Indices Monthly L3 Global 0.05Deg CMG (MOD13C2) ( MOD13C2网址 )，由"0.05Deg"可知该产品为经纬坐标系。

注意：暂时不推荐使用MATLAB处理MODIS的分幅正弦投影产品，即产品名称的空间分辨率为250 m、500 m、1 km等，建议用其它软件或GEE平台进行处理。

在直接读取MODIS产品前，我们可以先用 hdfinfo 来查看该数据的信息，注意括号内，最好加上'eos'，才能使读出来的信息容易辨读：

hdfinfo('.\Data\hdf_MOD13C2_NDVI\MOD13C2.A2021182.061.2021226060459.hdf', 'eos')

如图4所示，读取后的变量为struct类型，我们直接双击Grid，可以看到这里的Name字段 'MOD_Grid_monthly_CMG_VI' ，这是我们第一个要记录的东西。紧接着双击DataFields，如图5所示：

这里的Name字段，是我们要记录的第2个东西。例如我想读取的是MOD13C2的EVI数据，那我就得记录好'CMG 0.05 Deg Monthly NDVI'这个字段名。记录好后，我们就可以通过 hdfread 函数来读取出我们所要用到的数据：

data_hdf = hdfread('.\Data\hdf_MOD13C2_NDVI\MOD13C2.A2021182.061.2021226060459.hdf',...
        'MOD_Grid_monthly_CMG_VI', 'Fields', 'CMG 0.05 Deg Monthly NDVI');

括号内的 'MOD_Grid_monthly_CMG_VI' 即为Grid名， 'Fields' 为固定用法，千万别写成 'DataFields' 了， 'CMG 0.05 Deg Monthly NDVI' 便是DataField内的某一个数据集的名称。读取后的data_hdf为3600×7200的int16类型矩阵，我们简单使用下列函数看一下读取出的数据是否有误：

imagesc(data_hdf); colorbar;

如图6所示，2021年7月全球NDVI的空间分布便展示出来了，可能部分人不了解为什么这里有的值到了–2000，有的大于8000，这个就是后面《2 基于MATLAB的地理数据的预处理》的内容了。另外，通过上述命令可看出，这里的可视化并没有用到地理坐标的信息，因此imagesc给出的图仅仅是一幅图，而不是一幅带有地理信息的图像，那我们要怎么提取出该数据的地理坐标系呢？

根据图7官网给出的信息，我们可知该数据坐标系为地理经纬坐标系，覆盖范围为全球，即经度范围为–180°~180°，纬度范围为–90°~90°，分辨率为0.05°。因此我们可以推算出经纬网为：

lon_hdf = -180+0.025 : 0.05 : 180-0.025;
lat_hdf = 90-0.025 : -0.05 : -90+0.025;

注意，MATLAB中数据矩阵每一点的值代表的为地理栅格数据每个像元中心点的值，因此我们需要对一开始的边界，例如北纬90°，减去二分之一的空间分辨率即0.025°，来获得中心点的坐标。最终的经纬网可用meshgrid函数直接构建出来：

[lon_hdf, lat_hdf] = meshgrid(lon_hdf, lat_hdf);

如图8所示，data_hdf、lat_hdf和lon_hdf的矩阵大小完美契合，不放心的也可以双击打开lat_hdf和lon_hdf进行目视验证。关于如何对地理数据进行有地理信息的可视化这一部分的内容，将在《3 基于MATLAB的地理数据的可视化》这一章进行介绍。

(2) 基于MATLAB的SMAP数据的读取(HDF 5.0)

在 '\MATLAB_Course\Data\h5_SMAP\' 里，本文提供了2020年7月5日的全球土壤水分数据，所用产品为SMAP L3 Radiometer Global Daily 36 km EASE-Grid Soil Moisture, Version 8 (SPL3SMP) ( SPL3SMP网址 )，由"EASE-Grid"可知该产品为投影坐标系。

SMAP土壤水分产品的储存格式为HDF 5.0格式，具体体现在文件名的后缀为.h5，因此不能像MODIS一样使用hdfread等指令，而是用 h5info 、 h5disp 、 h5read 一类的函数。在读取SMAP数据前，先用 h5info 或 h5disp 查看下数据属性，本文以 h5info 为例：

h5info('.\Data\h5_SMAP\SMAP_L3_SM_P_20200705_R18290_002.h5')

结果图9所示，这里我们需要注意的' Groups' 的Name字段，以及 'Datasets' 的Name字段。接下来使用 h5read 分别读取SMAP降轨(AM)土壤水分、经纬度：

data_h5 = h5read('.\Data\h5_SMAP\SMAP_L3_SM_P_20200705_R18290_002.h5',...
        '/Soil_Moisture_Retrieval_Data_AM/soil_moisture');
lon_h5 = h5read('.\Data\h5_SMAP\SMAP_L3_SM_P_20200705_R18290_002.h5',...
        '/Soil_Moisture_Retrieval_Data_AM/longitude'); 
lat_h5 = h5read('.\Data\h5_SMAP\SMAP_L3_SM_P_20200705_R18290_002.h5',...
        '/Soil_Moisture_Retrieval_Data_AM/latitude');

括号里一定要注意，Groups名' Soil_Moisture_Retrieval_Data_AM' 前应有一斜杠 / ，后面分别接Datasets名 'soil_moisture' 、 'longitude' 和 'latitude' 。

读取后的data_h5为964×406的single类型矩阵，我们简单使用下列函数看一下读取出的数据是否有误：

imagesc(data_h5); colorbar;

可以发现，图10好像有点奇怪？其实这是EASE-Grid投影的问题。检查lat_h5和lon_h5我们可以发现，整个数据不是我们惯性思维的“上北下南左西右东”的方向，而是“左北右南上西下东”，在这里我们简单的对矩阵进行转置就可以使其变成“上北下南左西右东”。

data_h5 = data_h5'; lon_h5 = lon_h5'; lat_h5 = lat_h5';

如果你足够细心的话，你会发现EASE-Grid是投影坐标系，但SMAP产品提供的是latitude和longitude。这里的经纬度是根据投影坐标转换后的经纬度，所以并不是等间隔分布的(自行使用 diff 函数进行验证)，但这种“良心”的做法也使我们不需要去自己构建坐标系或者自己进行坐标系转换了。

1.4 网络通用数据格式(NetCDF)

在 '\MATLAB_Course\Data\nc_GPM\' 里，本文提供了2022年4月10日–4月15日的的全球降水数据，所用产品为GPM IMERG Early Precipitation L3 1 day 0.1 degree × 0.1 degree V06 (GPM_3IMERGDE) ( GPM网址 )，由"0.1 degree x 0.1 degree"可知该产品为经纬坐标系。

整体的读取方法还是和HDF格式差不多，主要的函数有 ncinfo 和 ncread ，这里不做详述，只列出代码：

data_nc = ncread('.\Data\nc_GPM\3B-DAY-E.MS.MRG.3IMERG.20220410-S000000-E235959.V06.nc4', 'precipitationCal');
lon_nc = ncread('.\Data\nc_GPM\3B-DAY-E.MS.MRG.3IMERG.20220410-S000000-E235959.V06.nc4', 'lon');
lat_nc = ncread('.\Data\nc_GPM\3B-DAY-E.MS.MRG.3IMERG.20220410-S000000-E235959.V06.nc4', 'lat');
lat_nc = flipud(lat_nc);