Как извлечь данные в файл netcdf на основе шейп-файла в R?

Question

У меня есть файл netcdf (глобальный домен), который я не знаю его проекционной информации, и я хочу извлечь данные (с lon и lat) на основе шейп-файла.

Цель состоит в том, чтобы найти данные для домена, представленного шейп-файлом (исходный файл netcdf содержит глобальные данные).Кроме того, окончательный формат данных должен быть фреймом данных, который содержит lon, lat и data.Я предполагаю, что функции extract и mask будут полезны.

Файлы 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/ 

Как мы видим, информация о проекции или crs для файла netcdf отсутствует.

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))

Проецировав оба на одну и ту же систему координат, мы можем применить шейп-файл в качестве маски к растровому кирпичу:

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