Incorrect length unit from read_ncdf

[admccurdy]

When importing a NetCDF file using read_ncdf the length unit of the imported object is incorrectly set as meter (EPSG: 9001) instead of kilometre (EPSG: 9036). If the file is read with either read_stars or terra::rast the length unit is correctly imported. I've only noticed the issue with files from the ORNL DAAC NetcdfSubset but I've only just started converting old code from raster to stars. Below is a toy example:

library(stars)

daymetr::download_daymet_ncss(
  location = c(41, -109, 36, -102),
  start = 1980,
  end = 1980,
  frequency = "annual",
  param = c("prcp"),
  path = getwd(),
  silent = F)  

dayMet_read_stars <- read_stars("prcp_annttl_1980_ncss.nc")
dayMet_read_ncdf <- read_ncdf("prcp_annttl_1980_ncss.nc")
dayMet_terra <- terra::rast("prcp_annttl_1980_ncss.nc")

st_crs(dayMet_read_stars)
st_crs(dayMet_read_ncdf)
terra::crs(dayMet_terra)

[edzer]

@dblodgett-usgs any ideas?

You probably want to read the precipitation, e.g. with

> read_stars("prcp_annttl_1980_ncss.nc", sub = 3)
prcp, 
stars object with 3 dimensions and 1 attribute
attribute(s):
            Min. 1st Qu. Median     Mean 3rd Qu.    Max.
prcp [mm] 132.75  301.59    362 411.8404  462.13 1839.85
dimension(s):
     from  to                  offset delta  refsys point values x/y
x       1 619                 -778.75     1 unnamed    NA   NULL [x]
y       1 569                  -119.5    -1 unnamed    NA   NULL [y]
time    1   1 1980-07-01 12:00:00 UTC    NA POSIXct    NA   NULL    

[admccurdy]

@edzer you are of course correct I would like to read the precipitation not lat :)

[dblodgett-usgs]

This is a tricky one. The data source has:

        short lambert_conformal_conic ;
                lambert_conformal_conic:latitude_of_projection_origin = 42.5 ;
                lambert_conformal_conic:false_easting = 0. ;
                lambert_conformal_conic:false_northing = 0. ;
                lambert_conformal_conic:standard_parallel = 25., 60. ;
                lambert_conformal_conic:semi_major_axis = 6378137. ;
                lambert_conformal_conic:inverse_flattening = 298.257223563 ;
                lambert_conformal_conic:grid_mapping_name = "lambert_conformal_conic" ;
                lambert_conformal_conic:longitude_of_central_meridian = -100. ;

the cf convention defines semi_major_axis in meters. There is nothing to say that the projection units are km.

However, the dataset has:

        float x(x) ;
                x:units = "km" ;
                x:long_name = "x coordinate of projection" ;
                x:standard_name = "projection_x_coordinate" ;

and

        float prcp(time, y, x) ;
                prcp:_FillValue = -9999.f ;
                prcp:missing_value = -9999.f ;
                prcp:coordinates = "time y x " ;
                prcp:grid_mapping = "lambert_conformal_conic" ;
                prcp:cell_methods = "area: mean time: sum within days time: sum over days" ;
                prcp:units = "mm" ;
                prcp:long_name = "annual total precipitation" ;

I can try to use the units of the x/y coordinate variables to define the units of the projection, but I'm afraid this will be pretty fragile. PR incoming.

[dblodgett-usgs]

With #534 I see:

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE

daymetr::download_daymet_ncss(
    location = c(41, -109, 36, -102),
    start = 1980,
    end = 1980,
    frequency = "annual",
    param = c("prcp"),
    path = getwd(),
    silent = T)  

dayMet_read_ncdf <- read_ncdf("prcp_annttl_1980_ncss.nc")
#> no 'var' specified, using prcp
#> other available variables:
#>  lat, y, x, lambert_conformal_conic, lon, time
#> Will return stars object with 352211 cells.
#> Warning in getGeoDatum(gm): Didn't find a longitude of prime meridian for datum,
#> assuming 0.

sf::st_crs(dayMet_read_ncdf)$ud_unit
#> 1 [km]

^{Created on 2022-05-04 by the reprex package (v2.0.1)}

[admccurdy]

Thanks @dblodgett-usgs and @edzer installed the PR and it works on my end, hopefully it's helpful for some additional use cases and isn't super specific to this dataset/method of access.