0

显然,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)

在此处输入图像描述

4

1 回答 1

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=立体 +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 的建议。

于 2020-05-11T04:24:25.340 回答