Forwarded from Ole Nielson on geonode-dev
Dear all,
We are developing an application for risk modelling based on spatial
data hosted in GeoServer and accessed through REST and OGC Services.
We came across what at first looked like precision issues in raster
resolutions when resolutions of downloaded data was different to what
was requested. This causes problems not just with our numerical
calculations but also for the integrity checks we have implemented
(e.g. that two layers have the same geo transform up to machine
precision).
However, a closer look reveals that it is probably more a matter of
what GeoServer does if asked for an impossible combination of bounding
box and grid resolution. The little study below would indicate that
the bounding box takes precedence and that resolution will be adjusted
from the requested resolution to something that is consistent.
Can someone please let us know if this assumption is correct and also
how to robustly compute consistent combinations of bounding
box/resolution so that the response from GeoServer is guaranteed to be
as requested?
Alternatively, what are your thoughts on changing the algorithm so
that the bounding box is changed to follow the grid point resolution
(I believe this is what is done in tools like Arc)?
Cheers and thanks
Ole
PS - I have attached a little Python program to assist with this study
--
Downloading using native parameters
--
Before uploading the file, gdalinfo reports
Origin = (122.5775001,2.0225000)
Pixel Size = (0.025,-0.025)
Corner Coordinates:
Upper Left ( 122.5775000, 2.0225000) (122d34'39.00"E, 2d 1'21.00"N)
Lower Left ( 122.5775000, -2.0025000) (122d34'39.00"E, 2d 0'9.00"S)
Upper Right ( 126.6025000, 2.0225000) (126d36'9.00"E, 2d 1'21.00"N)
Lower Right ( 126.6025000, -2.0025000) (126d36'9.00"E, 2d 0'9.00"S)
Center ( 124.590, 0.010) (124d35'24.00"E, 0d 0'36.00"N)
Downloading with REST using the exact same parameters
curl
"http://localhost:8001/geoserver-geonode-dev/ows?version=1.0.0&service=wcs&request=getcoverage&format=GeoTIFF&store=false&coverage=shakemap_20110505155015&crs=EPSG:4326&bbox=122.5775,-2.0025,126.6025,2.0225&resx=0.0250&resy=0.0250";
> dnld.tif
gdalinfo reports the exact same numbers as expected.
--
Downloading with adjusted bounding box
--
However, if I change the bounding box slightly (western boundary from
122.5775 to 122.58) and download again with the revised command:
curl
"http://localhost:8001/geoserver-geonode-dev/ows?version=1.0.0&service=wcs&request=getcoverage&format=GeoTIFF&store=false&coverage=shakemap_20110505155015&crs=EPSG:4326&bbox=122.58,-2.0025,126.6025,2.0225&resx=0.0250&resy=0.0250";
> dnld.tif
gdalinfo reports resolution values adjusted for the change in bounding box:
Origin = (122.578,2.0225000)
Pixel Size = (0.024984472049690,-0.025)
Corner Coordinates:
Upper Left ( 122.580, 2.0225000) (122d34'48.00"E, 2d 1'21.00"N)
Lower Left ( 122.580, -2.0025000) (122d34'48.00"E, 2d 0'9.00"S)
Upper Right ( 126.6025000, 2.0225000) (126d36'9.00"E, 2d 1'21.00"N)
Lower Right ( 126.6025000, -2.0025000) (126d36'9.00"E, 2d 0'9.00"S)
Center ( 124.5912500, 0.010) (124d35'28.50"E, 0d 0'36.00"N)
So instead of the requested 0.025 we get 0.02498447204969
-
Downloading with adjusted resolution
-
Similarly, and this is what triggered this question, if I change the
resolution slightly (from 0.025 to 0.03) but keep the original
bounding box and download with the command
curl
"http://localhost:8001/geoserver-geonode-dev/ows?version=1.0.0&service=wcs&request=getcoverage&format=GeoTIFF&store=false&coverage=shakemap_20110505155015&crs=EPSG:4326&bbox=122.5775,-2.0025,126.6025,2.0225&resx=0.0300&resy=0.0300";
> dnld.tif
gdal info reports
Origin = (122.5775001,2.0225000)
Pixel Size = (0.030037313432836,-0.030037313432836)
Corner Coordinates:
Upper Left ( 122.5775000, 2.0225000) (122d34'39.00"E, 2d 1'21.00"N)
Lower Left ( 122.5775000, -2.0025000) (122d34'39.00"E, 2d 0'9.00"S)
Upper Right ( 126.6025000, 2.0225000) (126d36'9.00"E, 2d 1'21.00"N)
Lower Right ( 126.6025000, -2.0025000) (126d36'9.00"E, 2d 0'9.00"S)
Center ( 124.590, 0.010) (124d35'24.00"E, 0d 0'36.00"N)
So instead of the requested 0.03 we get 0.030037313432836
"""Script to investigate relationship between specified bounding box and resolution in geoserver downloads via the REST interface.
"""
import os
bounding_box = [122.580, -2.0025
