如何根据R中的shapefile提取netcdf文件中的数据？

Question

我有一个 netcdf 文件（全局域），我不知道它的投影信息，我想根据 shapefile 提取数据（带有 lon 和 lat）。

目标是找到由 shapefile 表示的域的数据（原始 netcdf 包含全局数据）。此外，最终的数据格式应该是一个数据帧，其中包含 lon、lat 和 data。我认为extractandmask函数会很有用。

netcdf 和 shapefile 可以在https://www.dropbox.com/s/ju96eitzzd0xxg8/precip.2000.nc?dl=0和https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip? dl=0。谢谢你的帮助。

library(rgdal)
library(ncdf4)
library(raster)

shp = readOGR("gpr_000b11a_e.shp")
pre1 = nc_open("precip.2000.nc")
pre1

File precip.2000.nc (NC_FORMAT_NETCDF4_CLASSIC):
 1 variables (excluding dimension variables):
    float precip[lon,lat,time]   
        missing_value: -9.96920996838687e+36
        var_desc: Precipitation
        level_desc: Surface
        statistic: Total
        parent_stat: Other
        long_name: Daily total of precipitation
        cell_methods: time: sum
        avg_period: 0000-00-01 00:00:00
        actual_range: 0
         actual_range: 482.555358886719
        units: mm
        valid_range: 0
         valid_range: 1000
        dataset: CPC Global Precipitation

 3 dimensions:
    lat  Size:360
        actual_range: 89.75
         actual_range: -89.75
        long_name: Latitude
        units: degrees_north
        axis: Y
        standard_name: latitude
        coordinate_defines: center
    lon  Size:720
        long_name: Longitude
        units: degrees_east
        axis: X
        standard_name: longitude
        actual_range: 0.25
         actual_range: 359.75
        coordinate_defines: center
    time  Size:366   *** is unlimited ***
        long_name: Time
        axis: T
        standard_name: time
        coordinate_defines: start
        actual_range: 876576
         actual_range: 885336
        delta_t: 0000-00-01 00:00:00
        avg_period: 0000-00-01 00:00:00
        units: hours since 1900-01-01 00:00:00

7 global attributes:
    Conventions: CF-1.0
    version: V1.0
    history: created 9/2016 by CAS NOAA/ESRL PSD
    title: CPC GLOBAL PRCP V1.0
    References: https://www.esrl.noaa.gov/psd/data/gridded/data.cpc.globalprecip.html
    dataset_title: CPC GLOBAL PRCP V1.0
    Source: ftp://ftp.cpc.ncep.noaa.gov/precip/CPC_UNI_PRCP/ 

如我们所见，没有关于 netcdf 文件的投影或 crs 的信息。

Answer 1

您需要将 NetCDF 文件作为 raster* 对象打开，才能在 R 中将其作为光栅处理。为此，请使用brick或stack，而不是nc_open

pre1.brick = brick("precip.2000.nc")

您会注意到该文件使用气候科学中的常规约定，以 0 到 360 度范围内的值给出经度：

extent(pre1.brick)
# class       : Extent 
# xmin        : 0 
# xmax        : 360 
# ymin        : -90 
# ymax        : 90 

我们可以转换为常规的 -​​180 到 180 经度 usint rotate

pre1.brick = rotate(pre1.brick)

现在让我们看看我们的光栅和多边形文件的投影是什么：

crs(shp)
# +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
crs(pre1.brick)
# +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

您可以看到两者都使用经纬度坐标，但在基准面和椭圆方面有所不同。通常建议为了精度和速度，在可能的情况下，投影矢量数据而不是栅格，以使它们都在同一坐标系中：

shp = spTransform(shp, crs(pre1.brick))

将两者投影到同一个坐标系后，我们可以将 shapefile 作为遮罩应用到光栅砖上：

pre1.mask = mask(pre1.brick, shp)

并转换为 data.frame。这里我只演示第一层。如果您愿意，您可以通过[[1]]在以下行中删除来一次完成所有图层

pre1.df = as.data.frame(pre1.mask[[1]], xy=TRUE)

如果需要，您可以删除具有 NA 的行，只保留包含数据的掩码内的单元格：

pre1.df = pre1.df[complete.cases(pre1.df),]
head(pre1.df)
#            x     y X2000.01.01.00.00.00
# 10278 -81.25 82.75            0.2647110
# 10279 -80.75 82.75            0.2721323
# 10280 -80.25 82.75            0.2797630
# 10281 -79.75 82.75            0.2875668
# 10282 -79.25 82.75            0.2925712
# 10283 -78.75 82.75            0.2995636