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

Question

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

objective是为shapefile表示的域查找数据（原始netcdf包含全局数据）。此外，最终的数据格式应该是一个数据框，其中包含经度、纬度和数据。我假设 extract 和 mask 函数会有用。

netcdf 和 shapefile 可以在 https://www.dropbox.com/s/ju96eitzzd0xxg8/precip.2000.nc?dl=0 and 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 文件作为栅格*对象打开，以便在 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

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

How to extract the data in a netcdf file based on a shapefile in R?

r

extract

mask

spatial