How to reproject tiled data?

[khufkens]

This example should remove all doubt on how to reproject tiled Daymet data. The example downloads a single tile around Cape Cod, a distinct geographic feature. I plot the data in the corrected projection, as well as reprojected to EPSG 4326.

In short, the correct CRS projection of the tiled data is:
+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs

unlike what is reported in the tile meta-data!

library(daymetr)
library(raster)
library(sp)
library(rnaturalearth)

# download data
download_daymet_tiles(tiles = 11755, param = "prcp", start = 1980, end = 1980)

# read raster
r <- raster::raster(file.path(tempdir(),"prcp_1980_11755.nc"))

# set the correct projection
raster::projection(r) <- "+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"

# reproject to lat long
r_ll <- raster::projectRaster(r, crs = "+init=epsg:4326")

# grab usa outlines from rnaturalearth
usa_ll <- ne_countries(scale = 50, country = "United States of America", returnclass = "sp")
usa <- spTransform(usa_ll, CRSobj = CRS("+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))

# plot simple overview graphs
par(mfrow = c(2,1))
plot(r, main = "LLC projection")
lines(usa)

plot(r_ll, main = "lat lon projection")
lines(usa_ll)

[knussear]

I'm having the same problem - and this definitely does not solve it.

[khufkens]

When querying data using NCSS subsets (with coordinates, not tiled data) you might find that bounding boxes do not correspond to your query. You pass lat / lon coordinates which aren't orthogonal in the LCC "space". Hence, your bounding box of your query will change.

You find visual examples here:
https://daymet.ornl.gov/web_services.html

[khufkens]

For NCSS downloads the units of the projection are in kilometers not meters, use:

# set the correct projection
raster::projection(r) <- "+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=km +no_defs"