Setting ob_tran CRS #52
Comments
|
Currently, we can read this file as a curvilinear grid, but when read it as a regular GDAL file with |
|
> r = readGDAL("hmr_big.nc")
hmr_big.nc has GDAL driver netCDF
and has 1361 rows and 1286 columns
Warning message:
In readGDAL("hmr_big.nc") : GeoTransform values not availableso in order to test this, please provide a netcdf file that at least gives a longlat grid, not just a matrix, when read through GDAL. Or give the key to this (offset, delta for x and y). As far as I understand, we shouldn't read these from the long and lat matrices, because they would contain the non-obtran transformed coordinates (ie the curvilinear grid). |
|
Ok, I'll test this out and find a proper example in a few days. |
|
Here you can find a time slice of an example netCDF file from one EURO-CORDEX model. The grid is a rotated one, defined as:
Here is another similar example. Can this help? |
|
What would be the proj4string for these datasets? |
|
In theory it should be the same as previously discussed in r-spatial/sf#651, so something along the lines of |
|
This uses the forward projection trick also used in the issue you refer to. library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
(r = read_stars("example_ob_tran.nc"))
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# example_ob_tran.nc
# Min. :247.0
# 1st Qu.:267.4
# Median :275.9
# Mean :274.5
# 3rd Qu.:282.6
# Max. :292.6
# dimension(s):
# from to offset delta refsys point values
# x 1 106 -28.43 0.44 NA NA NULL [x]
# y 1 103 21.89 -0.44 NA NA NULL [y]
# time 1 1 1979-01-16 12:00:00 NA POSIXct NA NULL First the path using the lon = read_stars('NETCDF:"example_ob_tran2.nc":lon')
lat = read_stars('NETCDF:"example_ob_tran2.nc":lat')
rr = st_as_stars(r, curvilinear = list(x = lon[[1]], y = lat[[1]]))
rr
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# example_ob_tran.nc
# Min. :247.0
# 1st Qu.:267.4
# Median :275.9
# Mean :274.5
# 3rd Qu.:282.6
# Max. :292.6
# dimension(s):
# from to offset delta refsys point
# x 1 106 NA NA +proj=longlat +datum=WGS8... NA
# y 1 103 NA NA +proj=longlat +datum=WGS8... NA
# time 1 1 1979-01-16 12:00:00 NA POSIXct NA
# values
# x [106x103] -44.1407 [°], ..., 64.404 [°] [x]
# y [106x103] 22.1994 [°], ..., 72.4199 [°] [y]
# time NULL
# curvilinear grid
plot(rr, axes = TRUE, border = NA)
obtran = "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329"
st_crs(r) = 4326
r.ll = st_transform(r, obtran)
st_crs(r.ll) = 4326
plot(r.ll, axes = TRUE, border = NA)
r.ll
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# example_ob_tran.nc
# Min. :247.0
# 1st Qu.:267.4
# Median :275.9
# Mean :274.5
# 3rd Qu.:282.6
# Max. :292.6
# dimension(s):
# from to offset delta refsys point
# x 1 106 NA NA +proj=longlat +datum=WGS8... NA
# y 1 103 NA NA +proj=longlat +datum=WGS8... NA
# time 1 1 1979-01-16 12:00:00 NA POSIXct NA
# values
# x [106x103] -44.1407, ..., 64.404 [x]
# y [106x103] 22.1994, ..., 72.42 [y]
# time NULL
# curvilinear grid
st_bbox(rr)
# Units: [°]
# xmin ymin xmax ymax
# -44.14069 22.19937 64.40398 72.41995
st_bbox(r.ll)
# xmin ymin xmax ymax
# -44.14070 22.19937 64.40399 72.41996 |
|
I guess that now after #65 the another approach would be: |
|
the obtran PROJ string does not come from the file, |
|
I am having issues with fn = 'historical.nc'
fn %>% read_stars(sub = 'Q100', curvilinear=c('lon', 'lat'))
#Q100, trying to read file: netCDF:NETCDF:"/home/clima-archive4-b/afantini/regcm_simulations/EURO-CORDEX/change/italy_change/mrro/faster_way/data_hadgem/yearmax/historical.nc":Q100:lon
#Error in CPL_read_gdal(as.character(x), as.character(options), as.character(driver), :
file not found
fn %>% read_ncdf(var = 'Q100', curvilinear=c('lon', 'lat'))
#Error: length(d) == length(dim(x[[1]])) is not TRUE
lon = fn %>% read_ncdf(var = 'lon', quiet = TRUE)
lat = fn %>% read_ncdf(var = 'lat', quiet = TRUE)
ll = setNames(c(lon, lat), c("x", "y"))
fn %>% read_ncdf(var = 'Q100', quiet = TRUE) %>% st_as_stars(curvilinear = ll)
#Error: length(d) == length(dim(x[[1]])) is not TRUECan you reproduce this? |
fn = 'historical.nc'
lon = fn %>% read_ncdf(var = 'lon', quiet = TRUE)
lat = fn %>% read_ncdf(var = 'lat', quiet = TRUE)
lo = st_set_dimensions(lon, names = c("x", "y"))
la = st_set_dimensions(lat, names = c("x", "y"))
ll = setNames(c(lo, la), c("x", "y"))
x = fn %>% read_stars(sub = 'Q100', quiet = TRUE) %>% st_as_stars(curvilinear = ll)Any suggestions how we can make this easier? And, any idea what went wrong? |
|
Shouldn't we be able to read this by just doing |
|
If you read |
|
Thanks for sorting out! @mdsumner this may be related that in lon and lat both x and y dimensions have a negative delta (cellsize)? |
|
Is that the GDAL default? I remember it would always be upside down if the "default" transform was applied, where offsets both zero and deltas both 1. It was the wrong orientation for a simple png image. I'll explore |
|
First time I see a negative x delta! It is simply what the netcdf gdal driver communicates; I guess it comes from the ordering of the array in the netcdf file which is, afaict, not prescribed. |
fn = 'historical.nc'
fn %>% read_stars(sub = 'Q100', curvilinear = c('lon', 'lat')) %>% plotnow works. |
|
Mh, not for me with GDAL 2.2.2. I am getting: 2.2.2 |
|
What happens if you add |
|
Same error. |
|
This is after installing the branch with remotes::install_github("r-spatial/stars", "values")right? |
|
Dang! Sorry, I missed that. Yep, can confirm it's working if using that branch. I believe we can close this? |
|
Why don't we get the full time series in the stars object? |
|
Ah, sorry, that only concerns the |
|
Ah, see the damned 360-day calendar issue right there. |
|
I noticed - time for its own temporal reference system. Did you mention there was an R package to convert back and forth between calendar time and 360-day/year model time? |
|
@edzer where do you see a negative x? I see default geotransform on time, and typical +x/-y on the rest in gdalinfo terms Driver: netCDF/Network Common Data Format
Files: historical.nc
Size is 94, 114
Coordinate System is `'
Origin = (-9.509999923167690,-2.200000113090583)
Pixel Size = (0.109999998923271,-0.110000002700671)But, I'd assume to ignore the geotransform in curvilinear case anyway - right? (GDAL may or may not override the default, depending on vicarious heuristics). Here using raster (divorced from how I tackle this in ocean model data - ignoring the extent/crs, but analogous) and using the quadmesh grid/base plot: f <- "historical.nc"
r <- raster::raster(f, varname = "mrro")
coords <- raster::brick(raster::raster(f, varname = "lon"),
raster::raster(f, varname = "lat"))
library(quadmesh)
mesh_plot(r, coords = coords)
maps::map(add = TRUE)
|
Yup, |
|
@mdsumner I'm sorry, I mixed up the offset and delta colums! I tried to recreate the figure you made but failed and got this: |
|
Oh! Yes, so easy to mix up - worldfile order, GDAL order, raster extent order ... Yes both raster/quadmesh are CRAN versions. Hmm, I'm seeing that error in the netcdf pathway - so I'll investigate that(!): x <- read_ncdf(fn, curvilinear = c("lon", "lat"))
x
#Error in dim.stars(x) : length(d) == length(dim(x[[1]])) is not TRUE |
|
First shot at PCICt support: library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.7.0, GDAL 2.3.2, PROJ 5.2.0
fn = 'historical.nc'
(xx = read_stars(fn, sub = 'mrro', curvilinear = c("lon", "lat")))
# lon,
# lat,
# mrro,
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# historical.nc":mrro [kg/m^2/s]
# Min. :0
# 1st Qu.:0
# Median :0
# Mean :0
# 3rd Qu.:0
# Max. :0
# NA's :173400
# dimension(s):
# from to offset delta refsys point
# x 1 94 NA NA +proj=longlat +datum=WGS8... NA
# y 1 114 NA NA +proj=longlat +datum=WGS8... NA
# time 1 30 1976-07-01 360 days PCICt NA
# values
# x [94x114] 3.91253 [°],...,19.1685 [°] [x]
# y [94x114] 35.3469 [°],...,48.495 [°] [y]
# time NULL
# curvilinear grid
plot(xx[,,,1:12])
# Warning messages:
# 1: In seq(from = x$offset + (x$from - 1 + where) * x$delta, by = x$delta, :
# Incompatible methods ("+.PCICt", "Ops.difftime") for "+"
# 2: In seq(from = x$offset + (x$from - 1 + where) * x$delta, by = x$delta, :
# Incompatible methods ("+.PCICt", "Ops.difftime") for "+"Look at the labels: |
|
Wonderful! Should probably continue this in #29 |
I have a
netcdffile (link) of which I know the proper CRS, but the CRS is not encoded in the file itself:+proj=ob_tran +o_proj=longlat +o_lat_p=29.0 +o_lon_p=0.0 +lon_0=15.0 +a=6367470 +b=6367470stars, however, refuses to use it:This is very similar to what happened in r-spatial/sf#651 and has to do, as far as I understand it, with limited support in GDAL for the transformation
ob_tran.In the case of
sf, a workaround usinglwgeomwas possible. Is there anything similar that can be done forstars?The text was updated successfully, but these errors were encountered: