如何将水文降雨分析数据(MPEAHPS)栅格投影成可用的格式？

Question

显然，NOAA和NWS对他们的一些降雨数据采用了非传统的投影方式，在投影成传统格式方面并没有为其他用户提供很多帮助。我在让美国部分地区的栅格叠加方面取得了一点成功，但仍不太正确。

我希望有人能帮我解读一下我缺失的东西，并纠正这个数据的投影。

你可以在这里找到这个数据的更多信息。https: /polyploid.netblog?p=216https:/water.weather.govprecipdownload.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 +units=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. 你需要考虑是否需要双线或最近邻插值，如果你需要特定的输出栅格分辨率和范围，我建议你提供一个模型栅格，以用于 to 争论。

经编辑，结合@Robert Hijmans的建议。