Faulty second parallel as read by raster library

[khufkens]

- This is `raster` library bug, not a daymetr or ORNL one!!! 
- In short, raster fails to read the CF-1.6 netcdf format properly. 
- Only GDAL >=2.3 will do so properly through readGDAL().

The default projection data in returned tiles seems to be faulty. Mainly the second parallel North is set to 45 where it should be 60. Below is a reproducible example of the issue.

First download an example tile with a recognizable outline for reference

# set libraries and lat / lon projection
library(daymetr)
library(raster)
library(maps)
lat_lon <-  CRS("+proj=longlat +datum=WGS84")

# download a test tile with a faulty projection parameter
if(!file.exists("~/tmax_1980_11754.nc")){
  download_daymet_tiles(tile = "11754", path = "~", param = "tmax")
}

Read in the data with default projection info.

s <- raster("~/tmax_1980_11754.nc")

# print the default projection
projection(s)

Reproject this data to lat / lon and plot together with a lat / lon representation of the tiles and the coastal outline around NYC. Note the mismatch!

s_ll <- projectRaster(s, crs = lat_lon)
plot(s_ll)
plot(tile_outlines, add = TRUE)
maps::map("world",
          add = TRUE)

Now overwrite the original project with one where the second parallel North is 60 degrees instead of 45 (as specified on the website), and generate the plot again. Notice that all things align nicely now.

projection(s) <- CRS("+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0\
                     +y_0=0 +lat_1=25 +ellps=WGS84 +lat_2=60")
s_ll <- projectRaster(s, crs = lat_lon)

plot(s_ll)
maps::map("world",
     add = TRUE)
plot(tile_outlines, add = TRUE)

[khufkens]

The same issue exists for NCSS downloads

# download NCSS data
download_daymet_ncss(location = c(45.1, -70.1 , 45, -70),
                     start = 1980,
                     end = 1980,
                     param ="tmax",
                     path = "~")

# read in data
s <- raster("~/tmax_daily_1980_ncss.nc")

# plot projection, second parallel is 45 NOT 60
projection(s)

[khufkens]

Re-assign the projection of any tile (NCSS or otherwise) to fix the issue locally.

projection(s) <- CRS("+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +lat_1=25 +ellps=WGS84 +lat_2=60")

[khufkens]

After checking this seems to be a READ error not a write error while downloading from ORNL. This is a mess as this will either be an error in GDAL or NETCDF libraries!!

[khufkens]

Using GDAL 2.3 and either converting to geotiff or reading directly using

s <- rgdal::readGDAL("~/tmax_1980_11754.nc")

The coordinates read correctly:

Coordinate Reference System (CRS) arguments: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

However, the raster library messes things up!

[khufkens]

raster issue... closed for now