Open mikejohnson51 opened 2 years ago
hmm, quite right thank you! I see this a problem on Windows, but ok on Linux (only checked cursorily).
having a closer look, thanks
Awesome thanks @mdsumner! If useful, I am using a Mac.
ok, ta - can you run this for me? I think it's a {vapour} problem and this doesn't work correctly on Windows so if it's the same on Mac that's helpful (it should have different values in the two list elements, and one for each variable sst, ice, anom, err
vapour::vapour_sds_names("/vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/201809/oisst-avhrr-v02r01.20180929.nc")
but on Windows I only see
$datasource
[1] "/vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/201809/oisst-avhrr-v02r01.20180929.nc"
$subdataset
[1] "/vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/201809/oisst-avhrr-v02r01.20180929.nc"
I get:
vapour::vapour_sds_names("/vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/201809/oisst-avhrr-v02r01.20180929.nc")
#> $datasource
#> [1] "/vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/201809/oisst-avhrr-v02r01.20180929.nc"
#>
#> $subdataset
#> [1] "/vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/201809/oisst-avhrr-v02r01.20180929.nc"
Created on 2021-11-04 by the reprex package (v2.0.1)
ok thanks, now I'll check if it works with the local files (it should otherwise work well for other raster files you have, even online sources - but often they need extra syntax as you can see). Let me know if you have data sources you want to make work, it's unfortunate that my example no longer works (!)
Following on your examples, here is the data I am trying to get at (RAINRATE): https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/
(Unfortunately there is a 2 day rolling archive so that DateTime will not be available for long)...
vapour::vapour_sds_names
does pick up the sub datasets from the remote and local file.
The biggest difference is that its interpreting the remote subdatasets as HDF, and the local as NETCDF. Maybe its helpful?
library(gdalio)
#> Warning: package 'gdalio' was built under R version 4.1.1
url = 'https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc'
local = '/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc'
terra::rast(local)
#> class : SpatRaster
#> dimensions : 3840, 4608, 8 (nrow, ncol, nlyr)
#> resolution : 1000, 1000 (x, y)
#> extent : -2303999, 2304001, -1920000, 1920000 (xmin, xmax, ymin, ymax)
#> coord. ref. : Lambert_Conformal_Conic
#> sources : nwm.t00z.short_range.forcing.f001.conus.nc:U2D
#> nwm.t00z.short_range.forcing.f001.conus.nc:V2D
#> nwm.t00z.short_range.forcing.f001.conus.nc:LWDOWN
#> ... and 5 more source(s)
#> varnames : U2D (10-m U-component of wind)
#> V2D (10-m V-component of wind)
#> LWDOWN (Surface downward long-wave radiation flux)
#> ...
#> names : U2D, V2D, LWDOWN, RAINRATE, T2D, Q2D, ...
#> unit : m s-1, m s-1, W m-2, mm s^-1, K, kg kg-1, ...
#> time : 27266460
#### gdalio
grid0 <- list(extent = c( -2303999.25, 2304000.75, -1920000.375, 1919999.625),
dimension = c(3840, 4608),
projection = "+proj=lcc +lon_0=-90 +lat_1=33 +lat_2=45")
gdalio::gdalio_set_default_grid(grid0)
## REMOTE
(rmt = vapour::vapour_sds_names(glue::glue("/vsicurl/{url}", url = url)))
#> $datasource
#> [1] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [2] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [3] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [4] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [5] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [6] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [7] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [8] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#>
#> $subdataset
#> [1] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://LWDOWN"
#> [2] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://PSFC"
#> [3] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://Q2D"
#> [4] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://RAINRATE"
#> [5] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://SWDOWN"
#> [6] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://T2D"
#> [7] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://U2D"
#> [8] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://V2D"
gdalio_data(rmt$subdataset[4])
#> Warning in warp_in_memory_gdal_cpp(x, source_WKT = source_wkt, target_WKT = projection, : no valid projection in source, specified output projection will have no effect
#> (only the extent and dimension will be applied natively)
#> use 'source_projection' to apply the correct details for this source
#> we have a problem at RasterIO
#> $Band1
#> numeric(0)
### LOCAL
(loc = vapour::vapour_sds_names(local))
#> $datasource
#> [1] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [2] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [3] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [4] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [5] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [6] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [7] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [8] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#>
#> $subdataset
#> [1] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":U2D"
#> [2] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":V2D"
#> [3] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":LWDOWN"
#> [4] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":RAINRATE"
#> [5] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":T2D"
#> [6] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":Q2D"
#> [7] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":PSFC"
#> [8] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":SWDOWN"
str(gdalio_data(loc$subdataset[4]))
#> List of 1
#> $ Band1: num [1:17694720] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
Created on 2021-11-04 by the reprex package (v2.0.1)
ah, you'll have to use vrt() to tell vapour what the projection of the source is - is it "+proj=lcc +lon_0=-90 +lat_1=33 +lat_2=45" as you've used there?
GDAL sometimes doesn't pick that up, otherwise HDF5/NETCDF is kind of interchangeable - both remote and local SDS listing is looking fine
(thanks again)
it's working for me on Linux, so you might use
(ri <- vapour::vapour_raster_info(rmt$subdataset[4]))
$geotransform
[1] -2303999 1000 0 1920000 0 -1000
$dimXY
[1] 4608 3840
$minmax
[1] NA NA
$tilesXY
[1] 922 768
$projection
[1] "PROJCS[\"Lambert_Conformal_Conic\",GEOGCS[\"Unknown datum based upon the Authalic Sphere\",DATUM[\"Not_specified_based_on_Authalic_Sphere\",SPHEROID[\"Sphere\",6370000,0],AUTHORITY[\"EPSG\",\"6035\"]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],PARAMETER[\"central_meridian\",-97],PARAMETER[\"standard_parallel_1\",30],PARAMETER[\"standard_parallel_2\",60],PARAMETER[\"latitude_of_origin\",40],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]"
$bands
[1] 1
$projstring
[1] "+proj=lcc +lat_0=40 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs"
$nodata_value
and then
gdalio_data(vrt(rmt$subdataset[4], projection = ri$projection))
but, if ri$projection is missing, then you can add it manually here - I will try to find out why it's different on different machines :)
does that make sense?
Thank you for taking the time on this. Your suggestion makes sense and fixes my thinking (I think). I might still be running into a machine error which I am passing along.
In essence, the projection is not coming through in vapour::vapour_raster_info
. When added manually, it throws an error about The transformation is already "north up"...
. Adding transformation_options = 'SRC_METHOD=NO_GEOTRANSFORM'
get past the error, but only NaNs are returned. The same result is not seem on the local file:
Thanks again!
library(gdalio)
#> Warning: package 'gdalio' was built under R version 4.1.1
url = 'https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc'
(rmt = vapour::vapour_sds_names(glue::glue("/vsicurl/{url}", url = url)))
#> $datasource
#> [1] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [2] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [3] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [4] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [5] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [6] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [7] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [8] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#>
#> $subdataset
#> [1] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://LWDOWN"
#> [2] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://PSFC"
#> [3] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://Q2D"
#> [4] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://RAINRATE"
#> [5] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://SWDOWN"
#> [6] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://T2D"
#> [7] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://U2D"
#> [8] "HDF5:\"/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc\"://V2D"
(ri <- vapour::vapour_raster_info(rmt$subdataset[4]))
#> $geotransform
#> [1] 0 1 0 0 0 1
#>
#> $dimXY
#> [1] 4608 3840
#>
#> $minmax
#> [1] NA NA
#>
#> $tilesXY
#> [1] 922 768
#>
#> $projection
#> [1] ""
#>
#> $bands
#> [1] 1
#>
#> $projstring
#> [1] ""
#>
#> $nodata_value
#> [1] -9999
#>
#> $overviews
#> integer(0)
#>
#> $filelist
#> [1] "/vsicurl/https://nomads.ncep.noaa.gov/pub/data/nccf/com/nwm/prod/nwm.20211104/forcing_short_range/nwm.t00z.short_range.forcing.f001.conus.nc"
#>
#> $datatype
#> [1] "Float32"
#>
#> $extent
#> [1] 0 4608 3840 0
#proj string
prjstring = "+proj=lcc +lat_0=40 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs"
rainrate = gdalio_data(vrt(rmt$subdataset[4], projection = prjstring))
#> we have a problem at RasterIO
rainrate = gdalio_data(vrt(rmt$subdataset[4], projection = prjstring),
transformation_options = 'SRC_METHOD=NO_GEOTRANSFORM')
unique(rainrate$Band1)
#> [1] NaN
(loc = vapour::vapour_sds_names('/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc'))
#> $datasource
#> [1] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [2] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [3] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [4] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [5] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [6] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [7] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#> [8] "/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc"
#>
#> $subdataset
#> [1] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":U2D"
#> [2] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":V2D"
#> [3] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":LWDOWN"
#> [4] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":RAINRATE"
#> [5] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":T2D"
#> [6] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":Q2D"
#> [7] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":PSFC"
#> [8] "NETCDF:\"/Users/mjohnson/Downloads/nwm.t00z.short_range.forcing.f001.conus.nc\":SWDOWN"
rainrate = gdalio_data(vrt(loc$subdataset[4], projection = prjstring))
unique(rainrate$Band1)
#> [1] NaN 0.000000e+00 1.466858e-04 6.626712e-05 3.472222e-05
#> [6] 2.838869e-05 5.846876e-06 2.699782e-04 9.418342e-05 1.078746e-03
#> [11] 6.888158e-04 8.327957e-04 1.170947e-04 3.500261e-05 1.580749e-05
#> [16] 2.131663e-04 1.191407e-04 1.177377e-04 5.773942e-07 2.406020e-06
#> [21] 4.833812e-05 2.820079e-04 2.951334e-04 2.069941e-04 2.824237e-05
#> [26] 4.533716e-05 1.750545e-05 8.059172e-07 1.498722e-04 2.913366e-04
#> [31] 1.980195e-03 1.839105e-03 1.615666e-04 8.180914e-05 2.777778e-06
#> [36] 5.555556e-06 4.008175e-06 1.009084e-03 3.070762e-06 7.435790e-04
#> [41] 3.851855e-06 2.210039e-06 1.525097e-04 7.729674e-05 3.269683e-06
#> [46] 8.771006e-06 1.354910e-05 2.247516e-06 2.616020e-04 2.546528e-04
Created on 2021-11-05 by the reprex package (v2.0.1)
oh, the extent is wrong - it's just giving the degenerate index case 0, ncol, 0, nrow based on that default geotransform (0 1 0 0 0 1)
these files are just too problematic for simple use like this (simple being we can assume sensible extent, projection metadata is available)
you can wrap this stuff up in VRT, but it's quite specific how that gets done
on Linux this is working fine for me, though - no entirely sure why it's not working as well on Windows/mac
Hi @mdsumner. Yes being able to assume anything about these files (and even that the oddities don't change). Thanks for all the help. I will keep watching here to see if you can find any magic, or report if I am able to sort anything out.
I got quite lost with this, seeing differences whether it goes through the HDF5 or NetCDF driver, which I wouldn't have thought would make much difference - I'll let you know if I can figure it out, but really it's out of scope of gdalio. A future GDAL will have a very convenent and compact syntax for VRT, and that would make it easy to fix here.
Awesome - and thanks again. Happy to close/remain open depending on your preference
Hi! Thank for this package (among others). I am trying to run the Readme.Rmd example but get a persistent error of
fail to open
when accessing the remote NetCDF:I have tried,
I am not sure if this is (A) an issue in
gdalio
; (B) an issue with my install of GDAL or maybe (C) the URL format?I am hopeful you might have some insights on it ?
Below is a reprex. Thanks again for all you do.
Created on 2021-11-04 by the reprex package (v2.0.1)
PS: In the Readme.Rmd I think line 77 should read:
gdalio_data("my big.tif")