Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

netCDF in projected x,y reports bounds in getcaps as CRS:84 #68

Open
scaddenp opened this issue Sep 17, 2020 · 25 comments
Open

netCDF in projected x,y reports bounds in getcaps as CRS:84 #68

scaddenp opened this issue Sep 17, 2020 · 25 comments

Comments

@scaddenp
Copy link

Have netCDF which is stored in a projected x,y coordinates and I am pretty it is set up to CF standard. However, the getCaps call returns layers details where latlngBoundingBox = BoundingBox and BoundingBox is reported with SRS of CRS:84.

WMS calls with the requested SRS being the same as the native SRS are rather slow and I am wondering whether the data reader is converting from native grid to CRS:84 and then WMS is converting to again to the requested SRS?
@guygriffiths
Copy link
Collaborator

The GetCapabilities <EX_GeographicBoundingBox> needs to be in CRS:84, and I think that's the only bounding box reported by GetCapabilties?

I don't think there will be any native -> CRS:84 -> native conversion, it should just read the native grid. If you can send a file/request URL combination that's giving you trouble I can have a look at it.

@scaddenp
Copy link
Author

getCaps, in both 1.1.1 and 1.3.0 report

as well as Ex_geographicBoundingBox (1.3.0) or LatLonBoundingBox (1.1.x)

I have a backend mapping system which supports multiple different webservices (eg see https://data.gns.cri.nz/tez). Building it I have worked with many different implementations of WMS. in all of these, I successfully used either a dedicated tag or the SRS attribute or BoundingBox to determine the "native" SRS of a WMS layer. It isnt a big deal because it is only required for CQL queries with spatial predicates in systems that require those in native SRS. It did make my wonder though as to whether double projection was occuring. I want to look at the data extraction code anyway, so I will follow up on that.

@guygriffiths
Copy link
Collaborator

Ah, yes I see the issue, it's not differentiating between EX_GeographicBoundingBox and BoundingBox. Can you send me a data file which is giving you issues? I'm having trouble finding something in a different CRS (everything non-CRS:84 I have also includes 2d lat/lon, so it shows up as CRS:84).

@scaddenp
Copy link
Author

I have sent you link to file via email. Doesnt cross the antemedian. I was under impression all CF-compliant netCDF needed lat/lon? The header for this one is:
netcdf hpm_outputs {
dimensions:
time = 76 ;
lay = 8 ;
x = 501 ;
y = 302 ;
realn = 10 ;
variables:
int time(time) ;
time:standard_name = "time" ;
time:axis = "T" ;
time:long_name = "time" ;
time:units = "days since 1940-07-01 00:00:00" ;
char transverse_mercator ;
transverse_mercator:grid_mapping_name = "transverse_mercator" ;
transverse_mercator:longitude_of_central_meridian = 173. ;
transverse_mercator:false_easting = 1600000. ;
transverse_mercator:false_northing = 10000000. ;
transverse_mercator:latitude_of_projection_origin = 0. ;
transverse_mercator:scale_factor_at_central_meridian = 0.9996 ;
transverse_mercator:long_name = "CRS definition" ;
transverse_mercator:longitude_of_prime_meridian = 0. ;
transverse_mercator:semi_major_axis = 6378137. ;
transverse_mercator:semi_minor_axis = 6356752.31414036 ;
transverse_mercator:reference_ellipsoid_name = "GRS 1980" ;
transverse_mercator:prime_meridian_name = "Greenwich" ;
transverse_mercator:geographic_crs_name = "NZGD2000" ;
transverse_mercator:horizontal_datum_name = "New Zealand Geodetic Datum 2000" ;
transverse_mercator:projected_crs_name = "NZGD2000 / New Zealand Transverse Mercator 2000" ;
transverse_mercator:inverse_flattening = 298.257222101 ;
transverse_mercator:spatial_ref = "PROJCRS["NZGD2000 / New Zealand Transverse Mercator 2000",BASEGEOGCRS["NZGD2000",DATUM["New Zealand Geodetic Datum 2000",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4167]],CONVERSION["New Zealand Transverse Mercator 2000",METHOD["Transverse Mercator", ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]],PARAMETER["Longitude of natural origin",173,ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",1600000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",10000000,LENGTHUNIT["metre",1],ID["EPSG",8807]]], CS[Cartesian,2], AXIS["northing (N)",north,ORDER[1],LENGTHUNIT["metre",1]], AXIS["easting (E)",east,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["unknown"],AREA["New Zealand - onshore"],BBOX[-47.33,166.37,-34.1,178.63]],ID["EPSG",2193]]" ;
double lon(y, x) ;
lon:units = "degrees_east" ;
lon:long_name = "longitude coordinate" ;
lon:standard_name = "longitude" ;
double lat(y, x) ;
lat:units = "degrees_north" ;
lat:long_name = "latitude coordinate" ;
lat:standard_name = "latitude" ;
int lay(lay) ;
lay:standard_name = "model_level_number" ;
lay:long_name = "Model layer" ;
lay:positive = "down" ;
lay:axis = "Z" ;
lay:units = "1" ;
double x(x) ;
x:standard_name = "projection_x_coordinate" ;
x:long_name = "x coordinate of projection" ;
x:axis = "X" ;
x:units = "m" ;
double y(y) ;
y:standard_name = "projection_y_coordinate" ;
y:long_name = "y coordinate of projection" ;
y:axis = "Y" ;
y:units = "m" ;
int realn(realn) ;
realn:long_name = "Realisation" ;
realn:standard_name = "realization" ;
realn:units = "1" ;
float head(realn, time, lay, y, x) ;
head:_FillValue = -1.e+30f ;
head:grid_mapping = "transverse_mercator" ;
head:long_name = "Simulated GW Head" ;
head:standard_name = "(no standard name)" ;
head:units = "m" ;
head:coordinates = "lat lon" ;
float s_flow(realn, time, y, x) ;
s_flow:_FillValue = -1.e+30f ;
s_flow:grid_mapping = "transverse_mercator" ;
s_flow:long_name = "Simulated Stream Flow" ;
s_flow:standard_name = "(no standard name)" ;
s_flow:units = "m3/d" ;
s_flow:coordinates = "lat lon" ;
float s_flux(realn, time, y, x) ;
s_flux:_FillValue = -1.e+30f ;
s_flux:grid_mapping = "transverse_mercator" ;
s_flux:long_name = "Simulated Stream Flux to GW" ;
s_flux:standard_name = "(no standard name)" ;
s_flux:units = "m3/d" ;
s_flux:coordinates = "lat lon" ;

// global attributes:
:title = "hpm model outputs" ;
:description = "ensemble simulated outputs" ;
:file_creation_time = "2020-08-10 17:03:38.054010" ;
:Conventions = "CF-1.6" ;
:institution = "GNS Science" ;
:source = "PEST++/pyEMU" ;
}

@guygriffiths
Copy link
Collaborator

I am wondering whether the data reader is converting from native grid to CRS:84 and then WMS is converting to again to the requested SRS?

I've looked into it and yes, this is indeed what's happening. Transformed grids are also reporting their native CRS as CRS:84, which is what is appearing in the capabilities document.

The underlying issue here is that the NetCDF-Java library we are using gives us a co-ordinate system with a Projection object which we can use to translate between CRS:84 and the native projection of the data (and vice versa). But what it doesn't give us is a nicely encapsulated description of the native projection, nor an EPSG identifier. Without those things, it's not possible to properly report the native projection, nor determine whether co-ordinates are in the same reference system.

So for your example data, what ncWMS is essentially doing is:
Position in EPSG:8805 (from request URL) -> Position in CRS:84 -> Position in unidentified CRS (which happens to be EPSG:8805)

So it's not currently possible to fix this. I'm going to leave the issue open though, since the NetCDF-Java library is currently undergoing a fairly big revision, and there's a chance that this sort of thing may be supported in the not too distant future.

@scaddenp
Copy link
Author

scaddenp commented Sep 22, 2020 •
edited
Loading

Ok. Good to know. Perhaps I will raise issue on that library. I notice the dependency is 4.5.1 whereas netcdf-java is currently at 5.3.3. Maybe I should verify the problem on latest library first?

@ghansham
Copy link

ghansham commented Sep 22, 2020 via email

@scaddenp
Copy link
Author

Ghansham, thanks for that perspective. I dont think it matters whether WKT or EPSG (rather heavily used for many things) as Apache SIS (and other systems) can do either but the core issue for me, is that if data is stored in a projected coordinate system, then can it be extracted directly in the same coordinate system without code doing unnecessary projection/unprojections. I havent yet had time to delve into the reader code to see what goes on.

@ghansham
Copy link

ghansham commented Sep 23, 2020 via email

@guygriffiths
Copy link
Collaborator

Ok. Good to know. Perhaps I will raise issue on that library. I notice the dependency is 4.5.1 whereas netcdf-java is currently at 5.3.3. Maybe I should verify the problem on latest library first?

I'm going to do a release this week that depends on 5.3.3, but we've been using the 5.x branch for years - where are you seeing the 4.5.1 dependency?

@scaddenp
Copy link
Author

scaddenp commented Sep 23, 2020 via email

@guygriffiths
Copy link
Collaborator

Ah, you've probably seen that we're using version 5.1.0 of the netcdf4 library. The file is netcdf4-5.1.0.jar which scans as 4.5.1 at first (and often second and third!) glance.

@ghansham
Copy link

ghansham commented Sep 23, 2020 via email

@scaddenp
Copy link
Author

Ouch. You are right of course. it would have been more obvious if I had checked the pom.

Let me try and get a further understanding of the problem. There are two part to my mind:
1/ reporting the SRS in WMS
2/ extracting the data without reprojecting from the netcdf

1/ is a "nice to have". The Prj2EPSG API or the code from opengeo could do that. A rewrite would be desirable however to work with the Apache SIS resources.

2/ is more important. Do I understand that if KNOW the native coordinates in the netCDF are the same as the requested SRS, then it would be possible to extract directly? If so, then this should be achievable. A projection object can be created from requested SRS and components of projection object compared on a "fast fail" basis.

@ghansham
Copy link

ghansham commented Sep 23, 2020 via email

@scaddenp
Copy link
Author

I am pretty sure data is read once but projecting a entire grid twice is a computational cost. To me, the ncWMS is pretty slow compared to other WMS sources (we make heavy use of Geoserver and ArcGIS Server). Part of that could be improved with tile caching, but eliminating unnecessary calculations should also help.

@ghansham
Copy link

ghansham commented Sep 24, 2020 via email

@scaddenp
Copy link
Author

Got a hint on how to do that switch? is that fiddling with CdmUtils?

@ghansham
Copy link

ghansham commented Sep 24, 2020 via email

@scaddenp
Copy link
Author

Thanks. I will set up a timing test.

@guygriffiths
Copy link
Collaborator

Do I understand that if KNOW the native coordinates in the netCDF are the same as the requested SRS, then it would be possible to extract directly? If so, then this should be achievable. A projection object can be created from requested SRS and components of projection object compared on a "fast fail" basis.

Yes, essentially this behaviour is encapsulated in uk.ac.rdg.resc.edal.grid.cdm.CdmTransformedGrid#findIndexOf. It would need some modification, but if we had a method to compare SRS code with the Projection object and say whether they're equivalent, it'd be a straightforward change.

@scaddenp
Copy link
Author

Just getting back to this, and tried the change from scanline to bounding box. It is faster but not dramatically so. First time run (netCDF not in cache) it was 2.1sec per tile versus 2.4s for total load time of about 7.5s versus 8.4. Second run (netCDF in cache) it took total tile load time of 1.8sec versus 2.45. I suspect changing client to single tile load instead of tiles wms would be faster, but the antemeridian is likely a show stopper. I think a tile cache might be a better solution.

@ghansham
Copy link

ghansham commented Oct 22, 2020 via email

@scaddenp
Copy link
Author

Possibly - I admit that I havent peered closely at the file - I just used it for testing as I knew it was slow - and crossed the antemeridian. It is 600MB and can be found here https://share.gns.cri.nz/SPCKRK82KYT0/NZ_regional_BA_onshore_FAA_offshore_15as.nc.html

I could repeat tests with a smaller file, but I was looking for something that would highlight the difference.

@ghansham
Copy link

ghansham commented Oct 22, 2020 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@scaddenp @guygriffiths @ghansham