一、引言

在今天越来越严重的气候变化条件下，水文覆盖成为了越来越多研究者重视的话题。水文覆盖指的是地表或植被表面被水覆盖的面积，包括河流、洼地、湖泊、蓄水池等。它反应了一个地区的水资源分布、水域利用等情况，对于水资源管理和自然环境保护具有重要意义。

GRACE水文数据包括地表水蓄积（SWS）、土壤水蓄积（SSS）、总水蓄积（TWS）等变量，通常以每月为单位进行统计和融合，并以网格的形式提供各个区域的数据。

在这里，我们将通过使用 R 语言及其相关包对 GRACE 数据进行研究。具体来说，我们将使用 ncdf4 包读取 GRACE 的 .nc 数据，并进行数据的预处理和可视化分析。

二、数据获取

2.1 安装包和加载数据

# install.packages("ncdf4")
library(ncdf4)
ncdata <- nc_open("data.nc")
ncdata

2.2 基本数据信息解读

File D:\log\data\CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02.nc (NC_FORMAT_NETCDF4):

     2 variables (excluding dimension variables):
        float time_bounds[timebound,time]   (Contiguous storage)  
        float lwe_thickness[lon,lat,time]   (Contiguous storage)  
            grid_mapping: WGS84
            coordinates: time lat lon
            standard_name: Liquid_Water_Equivalent_Thickness
            long_name: Liquid_Water_Equivalent_Thickness
            Units: cm

     4 dimensions:
        timebound  Size:2 (no dimvar)
        time  Size:216 
            bounds: time_bounds
            calendar: gregorian
            axis: T
            standard_name: Time
            long_name: Time
            Units: days since 2002-01-01T00:00:00Z
        lon  Size:1440 
            bounds: lon_bounds
            valid_max: 359.875
            valid_min: 0.125
            axis: X
            standard_name: Longitude
            long_name: Longitude
            Units: degrees_east
        lat  Size:720 
            bounds: lat_bounds
            valid_max: 89.875
            valid_min: -89.875
            axis: Y
            standard_name: Latitude
            long_name: Latitude
            Units: degrees_north

    58 global attributes:
        Conventions: CF-1.6, ACDD-1.3, ISO 8601
        filename: netcdf/CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02.nc
        Metadata_Conventions: Unidata Dataset Discovery v1.0
        ......

从文件信息中，我们可以看出该文件包含四个变量：time（时间）、lon（经度）、lat（纬度）和 lwe_thickness（覆盖厚度）。

2.3 基础数据提取

lon <- ncvar_get(ncdata,'lon') # 经度
nrow(lon)
lat <- ncvar_get(ncdata,'lat') # 纬度
nrow(lat)

time <- ncvar_get(ncdata,'time') # 时间
time
lwe_thickness <- ncvar_get(ncdata,'lwe_thickness') # 覆盖厚度
class(lwe_thickness) # 判定数据类型

结果展示：

# 经度
[1] 1440
# 纬度
[1] 720

# 时间
[1]  107.0  129.5  227.5  258.0  288.5  319.0  349.5  380.5  410.0  439.5
 [11]  470.0  495.5  561.5  592.5  623.0  653.5  684.0  714.5  736.5  777.0
 [21]  805.5  836.0  866.5  897.0  927.5  958.5  989.0 1019.5 1050.0 1080.5
 [31] 1111.5 1141.0 1170.5 1201.0 1231.5 1262.0 1292.5 1323.5 1354.0 1384.5
 [41] 1415.0 1445.5 1476.5 1506.0 1535.5 1566.0 1596.5 1627.0 1657.5 1688.5
 [51] 1719.0 1749.5 1780.0 1810.5 1841.5 1870.5 1900.5 1931.0 1961.5 1992.0
 [61] 2022.5 2053.5 2084.0 2114.5 2145.0 2175.5 2206.5 2236.5 2266.5 2297.0
 [71] 2327.5 2358.0 2388.5 2419.5 2450.0 2480.5 2511.0 2541.5 2572.5 2602.0
 [81] 2631.5 2662.0 2692.5 2723.0 2753.5 2784.5 2815.0 2845.5 2876.0 2906.5
 [91] 2937.5 2967.0 2996.5 3027.0 3057.5 3088.0 3118.5 3149.5 3180.0 3210.5
[101] 3241.0 3269.5 3335.5 3361.5 3392.0 3422.5 3485.5 3514.5 3545.0 3575.5
[111] 3590.5 3648.0 3667.5 3697.5 3727.5 3746.0 3819.0 3849.5 3880.5 3908.5
[121] 3974.5 4002.5 4033.5 4062.0 4128.0 4153.5 4184.0 4214.5 4306.5 4337.0
[131] 4367.5 4391.0 4458.5 4488.0 4518.5 4546.0 4610.5 4641.0 4671.5 4702.0
[141] 4769.5 4793.0 4822.5 4853.0 4864.0 4943.5 4975.5 5004.5 5104.0 5128.5
[151] 5157.0 5188.5 5253.0 5280.0 5309.5 5346.5 5444.5 5471.5 5499.0 5568.5
[161] 5592.5 5611.0 5639.5 6010.0 6034.0 6147.5 6163.0 6193.5 6224.5 6254.0
[171] 6283.5 6314.0 6344.5 6375.0 6405.5 6436.5 6467.0 6497.5 6528.0 6558.5
[181] 6589.5 6619.5 6649.5 6680.0 6710.5 6741.0 6771.5 6802.5 6833.0 6863.5
[191] 6894.0 6924.5 6955.5 6985.0 7014.5 7045.0 7075.5 7106.0 7136.5 7167.5
[201] 7198.0 7228.5 7259.0 7289.5 7320.5 7350.0 7379.5 7410.0 7440.5 7471.0
[211] 7501.5 7532.5 7563.0 7593.5 7624.0 7654.5


# 覆盖厚度
[1] "array"

2.4 数据预处理

日期格式转换

# 这份测绘数据起始时间是2002.01.01，time的值表示过去的天数
time <- as.Date("2002-01-01") + time
time

结果展示：

[1] "2002-04-18" "2002-05-10" "2002-08-16" "2002-09-16" "2002-10-16"
  [6] "2002-11-16" "2002-12-16" "2003-01-16" "2003-02-15" "2003-03-16"
 [11] "2003-04-16" "2003-05-11" "2003-07-16" "2003-08-16" "2003-09-16"
 [16] "2003-10-16" "2003-11-16" "2003-12-16" "2004-01-07" "2004-02-17"
 [21] "2004-03-16" "2004-04-16" "2004-05-16" "2004-06-16" "2004-07-16"
 [26] "2004-08-16" "2004-09-16" "2004-10-16" "2004-11-16" "2004-12-16"
 [31] "2005-01-16" "2005-02-15" "2005-03-16" "2005-04-16" "2005-05-16"
 [36] "2005-06-16" "2005-07-16" "2005-08-16" "2005-09-16" "2005-10-16"
 [41] "2005-11-16" "2005-12-16" "2006-01-16" "2006-02-15" "2006-03-16"
 [46] "2006-04-16" "2006-05-16" "2006-06-16" "2006-07-16" "2006-08-16"
 [51] "2006-09-16" "2006-10-16" "2006-11-16" "2006-12-16" "2007-01-16"
 [56] "2007-02-14" "2007-03-16" "2007-04-16" "2007-05-16" "2007-06-16"
 [61] "2007-07-16" "2007-08-16" "2007-09-16" "2007-10-16" "2007-11-16"
 [66] "2007-12-16" "2008-01-16" "2008-02-15" "2008-03-16" "2008-04-16"
 [71] "2008-05-16" "2008-06-16" "2008-07-16" "2008-08-16" "2008-09-16"
 [76] "2008-10-16" "2008-11-16" "2008-12-16" "2009-01-16" "2009-02-15"
 [81] "2009-03-16" "2009-04-16" "2009-05-16" "2009-06-16" "2009-07-16"
 [86] "2009-08-16" "2009-09-16" "2009-10-16" "2009-11-16" "2009-12-16"
 [91] "2010-01-16" "2010-02-15" "2010-03-16" "2010-04-16" "2010-05-16"
 [96] "2010-06-16" "2010-07-16" "2010-08-16" "2010-09-16" "2010-10-16"
[101] "2010-11-16" "2010-12-14" "2011-02-18" "2011-03-16" "2011-04-16"
[106] "2011-05-16" "2011-07-18" "2011-08-16" "2011-09-16" "2011-10-16"
[111] "2011-10-31" "2011-12-28" "2012-01-16" "2012-02-15" "2012-03-16"
[116] "2012-04-04" "2012-06-16" "2012-07-16" "2012-08-16" "2012-09-13"
[121] "2012-11-18" "2012-12-16" "2013-01-16" "2013-02-14" "2013-04-21"
[126] "2013-05-16" "2013-06-16" "2013-07-16" "2013-10-16" "2013-11-16"
[131] "2013-12-16" "2014-01-09" "2014-03-17" "2014-04-16" "2014-05-16"
[136] "2014-06-13" "2014-08-16" "2014-09-16" "2014-10-16" "2014-11-16"
[141] "2015-01-22" "2015-02-15" "2015-03-16" "2015-04-16" "2015-04-27"
[146] "2015-07-15" "2015-08-16" "2015-09-14" "2015-12-23" "2016-01-16"
[151] "2016-02-14" "2016-03-16" "2016-05-20" "2016-06-16" "2016-07-15"
[156] "2016-08-21" "2016-11-27" "2016-12-24" "2017-01-21" "2017-03-31"
[161] "2017-04-24" "2017-05-13" "2017-06-10" "2018-06-16" "2018-07-10"
[166] "2018-10-31" "2018-11-16" "2018-12-16" "2019-01-16" "2019-02-15"
[171] "2019-03-16" "2019-04-16" "2019-05-16" "2019-06-16" "2019-07-16"
[176] "2019-08-16" "2019-09-16" "2019-10-16" "2019-11-16" "2019-12-16"
[181] "2020-01-16" "2020-02-15" "2020-03-16" "2020-04-16" "2020-05-16"
[186] "2020-06-16" "2020-07-16" "2020-08-16" "2020-09-16" "2020-10-16"
[191] "2020-11-16" "2020-12-16" "2021-01-16" "2021-02-15" "2021-03-16"
[196] "2021-04-16" "2021-05-16" "2021-06-16" "2021-07-16" "2021-08-16"
[201] "2021-09-16" "2021-10-16" "2021-11-16" "2021-12-16" "2022-01-16"
[206] "2022-02-15" "2022-03-16" "2022-04-16" "2022-05-16" "2022-06-16"
[211] "2022-07-16" "2022-08-16" "2022-09-16" "2022-10-16" "2022-11-16"
[216] "2022-12-16"

现在就是问我们能看懂的时间类型了，仔细看，这数据是每个月记录测绘一次结果，其中有不少缺失数据。下一篇我们将介绍如何补全这些测绘数据，欢迎关注和一起讨论。

lwe_thickness理解

lwe_thickness[lon,lat,time]数据是float类型，表示每个时间点的经度和纬度的覆盖厚度值。这里补充说明一下：一个数据矩阵代表一个时间节点的一片区域的水文覆盖厚度情况。了解这个，我们下一步就可以将栅格化数据转换成dataframe.

# 检查数组的维度
dim(lwe_thickness)

结果展示：

[1] 1440  720  216

是不是和lon、lat和time的长度是一致的。

扁平化数据

# 创建一个空的数据框，指定 time 为 Date 类型

for(i in 1:dim(lwe_thickness)[1]){
  for(j in 1:dim(lwe_thickness)[2]){
    for(k in 1:dim(lwe_thickness)[3]){
      # 添加一行数据
      new_row <- data.frame(time = time[k],lon = lon[i],lat = lat[j],lwe_thickness = lwe_thickness[i,j,k])
      df <- rbind(df, new_row)
    }
  }
}

library(data.table)

# 对数据进行重组和转置，以扁平化数据结构
lwe_thickness_flat <- as.data.table(lwe_thickness)

结果展示：

             V1  V2  V3     value
        1:    1   1   1 -2.468668
        2:    1   1   2 -1.569404
        3:    1   1   3 -3.044039
        4:    1   1   4 -3.201635
        5:    1   1   5 -1.785586
       ---                       
223948796: 1440 720 212  2.893622
223948797: 1440 720 213  5.812358
223948798: 1440 720 214  9.525175
223948799: 1440 720 215 10.497257
223948800: 1440 720 216  8.949245

数据量很恐怖呀，达到了惊人的2亿多条！

将V1、V2、V3分别替换成lon、lat和time的数据

# lon
lwe_thickness_flat$V1 <- lon[lwe_thickness_flat$V1]
# lat
lwe_thickness_flat$V2 <- lat[lwe_thickness_flat$V2]
# time
lwe_thickness_flat$V3 <- time[lwe_thickness_flat$V3]

修改列名

# 将数据框lwe_thickness_flat转换为数据表对象
lwe_thickness_flat<-as.data.table(lwe_thickness_flat)
# 设置列名
setnames(lwe_thickness_flat, c("V1","V2","V3","value"), c("lon","lat","time","lwe_thickness"))

由于数据量巨大，有2亿多条数据，无法直接修改列名。至此，数据已经提取完全。

存入csv

# 写入文件并忽略行名
fwrite(lwe_thickness_flat, file = "D:/log/modis.csv", row.names = FALSE)

转存成csv，免得以后分析的时候需要再次转化！

在这里插入图片描述

三、基础数据展示

抽样

由于数据量过于庞大，选取一个栅格的数据。

library(dplyr)

# 对数据进行日期筛选
sampled_data <- lwe_thickness_flat %>% 
                     filter(time == as.Date("2002-04-18"))

绘图


# 设定等值线的级别序列
level_seq <- seq(from = -90, to = 300, by = 0.25)

# 绘制水平等值线图，并使用上面设定的等值线级别
p <- ggplot(data = sampled_data, aes(x = lon, y = lat, z = lwe_thickness)) + 
    geom_contour(aes(colour = ..level..), breaks = level_seq) +
    scale_colour_gradient(low = "blue", high = "red")
# 将ggplot对象换成plotly对象
ggplotly(p)