load CDM CRS if they are available in string form #936
Conversation
|
Looks like I did something fishy here that leads CRS to return a bool |
|
When I test on my end I still get an using Rasters
import NCDatasets, Downloads
paths2granule = "https://hyp3-jhk-sandbox-contentbucket-wxvu5zy6v7d6.s3.us-west-2.amazonaws.com/100f0e79-d569-42b9-8a20-3c3384b34dee/S1B_IW_SLC__1SSH_20181227T194309_20181227T194320_014232_01A75C_12D4_X_S1A_IW_SLC__1SSH_20190102T194350_20190102T194401_025303_02CCA3_D38F_G0120V02_P099.nc"
local_path = joinpath(tempdir(), splitpath(paths2granule)[end])
isfile(local_path) || Downloads.download(paths2granule, local_path)
rs = RasterStack(local_path; name=[:v])
Rasters.crs(rs) |
|
So it looks like the RasterStack constructor doesn't follow this pipeline. julia> crs(Raster(local_path; name=:v))
ESRIWellKnownText with CRS mode: PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH],AUTHORITY["EPSG","3413"]] |
because WKT2 is guaranteed to be backwards compatible with WKT1
| @@ -339,6 +343,48 @@ function _cdmlookup( | |||
| end | |||
| return _cdmlookup(D, index, order, span, sampling, metadata, crs, mappedcrs) | |||
| end | |||
|
|
|||
| function _get_crs_by_hook_or_crook(ds::AbstractDataset, var, attr) | |||
There was a problem hiding this comment.
One problem with this is we are assuming which dimensions crs should attach to instead of reading the explicit link to x/y etc dimensions in the grid mapping.
I think we need to return something like ["lat", "lon"] => crs so we can match them properly when we make dimensions.
There was a problem hiding this comment.
explicit link to x/y etc dimensions in the grid mapping.
Isn't that in the actual x/y dims though? This is reading from the variable var primarily...
There was a problem hiding this comment.
Yeah but it still lists which dimensions of the variable the crs relates to, and we should probably keep that information even though nearly always its obvious which those should be.
There was a problem hiding this comment.
We don't have anywhere to store that in Projected though. And we already do that when choosing the dims for the variable right?
Or is that just following the order of definition?
That sounds too complicated and like too much of a rework to go in this PR though, if we just want CRS support. It can come later but I have a feeling that no one actually screws around with dimensions that much.
There was a problem hiding this comment.
Everything else about dimensions come just from the dimension variable and other variables (like bounds) that it explicitly links to.
This is quite different - its putting some metadata from another variable into the dimensions - so we need to make a link from variable to its dimensions - but some of those are the crs dimensions, some are other dimensions.
It is actually explicitly specified which dimensions the crs relates to, so instead of guessing (very educated guessing for sure) the link we should use that information.
There was a problem hiding this comment.
Yeah totally at the variable dimension level!
One reason is just to be defensive against unknown possibilities that follow the spec but are "weird" or unintuitive.
The current hack is marginally justifiable in that we can't read the metadata at all. But if we are reading it a bit, we should do it properly.
Unfortunately I think we need to read the CF standards and really get clear on how grid_mapping works ;)
There was a problem hiding this comment.
This is where we link dimensions to variables
https://github.com/rafaqz/Rasters.jl/blob/main/src/sources/commondatamodel.jl#L184-L189
CommonDataModel provides CDM.dimnames(var) that gives you the string names of the dimension variables connected to any variable.
There was a problem hiding this comment.
From the standard:
"An application can determine that x and y are the projection coordinates by recognizing the standard names projection_x_coordinate and projection_y_coordinate"
There was a problem hiding this comment.
This will link arbitrary dimension names to spatial x/y
(we should also use it to determine which are the X and Y dims...)
There was a problem hiding this comment.
reading this code again I finally get what you mean :D - we just need to make the dimension type lookup do something different from what it currently does (name pattern matching)
|
Should we get this in for now, if we can fix the issues with returning Bools, and then get the CF projection stuff out with the dimension linking whenever it's ready? CF also has different names depending on geodetic/projected CRS for dimensions, which is another pain... |
|
No... Because then it will be breaking to do the dimensions linking ;) |
|
Wouldn't that be breaking now as well? Or is it not breaking because we don't return anything now? |
|
Yeah maybe they're both breaking |
we still cannot load CF CRS natively, but most writers will also write the CRS in string form, which we can read. The key indicates the format.
Fix #935
Partially fix #31