Как спроецировать растр Hydrologi c Data Analysis Data (MPE / AHPS) в пригодный для использования формат?

Question

Очевидно, NOAA и NWS используют нетрадиционный прогноз для некоторых своих данных об осадках и не предлагают большой помощи с точки зрения проецирования их в традиционный формат для других пользователей. Я добился небольшого успеха в наложении растра для части Соединенных Штатов, но это все еще не совсем то.

Я надеюсь, что кто-нибудь поможет мне расшифровать то, что мне не хватает, и исправить проекцию этих данных.

Вы можете найти более подробную информацию об этих данных здесь: https://polyploid.net/blog/?p=216 https://water.weather.gov/precip/download.php

library(tidyverse)
library(raster)
library(rgdal)
library(sp)

setwd("C:/Users/MPE_Data/")

file_list <- list.files("201809")

grib0<-raster::brick("201809//ST4_2018091307_24h.nc", varname="APCP_SFC")[[1]]
grib0@crs
crs(grib0) <- "+proj=longlat +a=6371200 +b=6371200 +no_defs"
crs(grib0) <- "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-105 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs"


us_shp <- rgdal::readOGR("C:/Users/cb_2017_us_state_500k/US_clipped.shp")
shp <- rgdal::readOGR("C:/Users/nc_sc_counties_wgs1984.shp")

wgs<-"+proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs"
wgsraster <- projectRaster(grib0, crs=wgs)
plot(wgsraster)
shp <- spTransform(shp, CRS(wgs))
us_shp <- spTransform(us_shp, CRS(wgs))
plot(shp,add=TRUE)
plot(us_shp,add=TRUE)

Answer 1

Мне не удалось найти вашу точную карту, но вот пример с использованием последних данных об осадках. Вам не нужно назначать CRS, поскольку файл netCDF уже имеет связанный с ним CRS, вы можете просто projectRaster. Также на веб-сайте NOAA есть возможность загрузки в geoTIFF, которую я бы порекомендовал, если вам это удобнее.

require(raster)
require(ncdf4)
require(maptools)

data(wrld_simpl)
us_shp=wrld_simpl[which(wrld_simpl$NAME=="United States"),]

rs=raster::brick("./nws_precip_1day_20200509_netcdf/nws_precip_1day_20200509_conus.nc",varname="observation")[[1]]
rs@crs ##note already has a crs associated with it

+ proj = stere + lat_0 = 90 + lat_ts = 60 + lon_0 = -105 + x_0 = 0 + y_0 = 0 + a = 6371200 + b = 6371200 + единицы = m + no_defs

##assign the pixels with -10000 to NA.    
NAvalue(rs) = -10000

##reproject to longlat WGS84
rs=projectRaster(rs,crs=crs(us_shp))
plot(rs,col=rainbow(100))
lines(us_shp)
##note the data extends outside the bounds of country

##use mask to remove data that is not over the land area
rs=mask(rs,us_shp)
plot(rs,col=rainbow(100)
lines(us_shp)

введите описание изображения здесь

Обратите внимание, что максимальное значение rs изменилось с 7,8 на 7,0 из-за метода билинейной интерполяции, используемого в projectRaster. Вам нужно решить, требуется ли вам билинейная интерполяция или интерполяция ближайшего соседа, и если вам нужно указать c разрешение и размер выходного растра, я бы предложил предоставить растр модели для аргумента to.

Отредактировано с учетом предложения @Robert Hijmans.