How to project Hydrologic Rainfall Analysis Data (MPE/AHPS) raster to a usable format?

Question

Apparently NOAA and the NWS use a non-traditional projection for some of their rainfall data and don't offer a lot of help in terms of projecting it to a traditional format for other users. I've had a bit of success in getting the raster to overlay for part of the United States but it still isn't quite right.

I'm hoping someone can help me decipher what I am missing and correct the projection of this data.

You can find more information of this data here: 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

I couldn't find your exact map but here is an example using recent precipitation data. You don't need to assign a CRS as the netCDF file already has a CRS associated with it, you can simply projectRaster. Also the NOAA website has the option to download to geoTIFF which I would recommend if you are more comfortable with that.

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)

Note that the maximum value of rs changed from 7.8 to 7.0 due to the bilinear interpolation method used in projectRaster. You need to consider whether you require bilinear or nearest neighbour interpolation and if you need to be specific about the output raster resolution and extent I would suggest supplying a model raster for the to argument.

Edited to incorporate @Robert Hijmans' suggestion.