Opened 16 years ago
Last modified 8 years ago
#352 new defect
r.in.gdal region troubles in LL
Reported by: | neteler | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | 6.4.6 |
Component: | Raster | Version: | unspecified |
Keywords: | r.in.gdal | Cc: | Agustin.Lobo@… |
CPU: | Unspecified | Platform: | Unspecified |
Description
The GeoTIFF file http://aloboaleu.googlepages.com/npp.tif (3.3MB) is causing problems in the region settings when imported with r.in.gdal:
gdalinfo npp.tif Driver: GTiff/GeoTIFF Files: npp.tif Size is 1440, 572 Coordinate System is: GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.2572235630016, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","4326"]] Origin = (-180.000000000000000,85.000000000114397) Pixel Size = (0.250000000000200,-0.250000000000200) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left (-180.0000000, 85.0000000) (180d 0'0.00"W, 85d 0'0.00"N) Lower Left (-180.0000000, -58.0000000) (180d 0'0.00"W, 58d 0'0.00"S) Upper Right ( 180.0000000, 85.0000000) (180d 0'0.00"E, 85d 0'0.00"N) Lower Right ( 180.0000000, -58.0000000) (180d 0'0.00"E, 58d 0'0.00"S) Center ( 0.0000000, 13.5000000) ( 0d 0'0.00"E, 13d30'0.00"N) Band 1 Block=1440x1 Type=Float32, ColorInterp=Gray g.region n=90 s=-90 w=-180 e=180 res=0:05 -p projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 90N south: 90S west: 180W east: 180E nsres: 0:05 ewres: 0:05 rows: 2160 cols: 4320 cells: 9331200 r.in.gdal npp.tif out=npp Projection of input dataset and current location appear to match 100% <npp> created r.in.gdal complete. g.region rast=npp -p projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 85N south: 58S west: 180W east: 179:59:59.999999W nsres: 0:15 ewres: 0 rows: 572 cols: 1440 cells: 823680 g.region -p ERROR: region for current mapset is invalid line 9: <e-w resol: 0> run "g.region"
Apparently the global wrap-around recognition is failing for 180W/E.
Markus
Change History (6)
follow-up: 2 comment:1 by , 16 years ago
follow-up: 3 comment:2 by , 16 years ago
Cc: | added |
---|
Replying to glynn:
Replying to neteler:
gdalinfo npp.tif
Pixel Size = (0.250000000000200,-0.250000000000200)
It's that extra 0.0000000000002 that's causing the problems:
I see.
...
[And which will forever cause problems. We'll just have to work around them as we find them. Almost everything that you "know" about coordinate geometry is wrong for spherical geometry. E.g. just because x1 == x2, that doesn't mean that x1-x0 == x2-x0.]
The user took this map
http://user.uni-frankfurt.de/~grieser/downloads/NetPrimaryProduction/npp_CruP_All_fine_geo.zip
and passed it to R, and then redefined:
>b <- npp_geotiff >b@data$band1[b@data$band1 <0 ]<- 0
and then he wrote b as npp.tif using the QGIS plugin manageR.
(Source: http://lists.osgeo.org/pipermail/grass-user/2008-October/047211.html)
Maybe there is a precision issue before, somwhere in the R or QGIS part?
Markus
follow-up: 6 comment:3 by , 16 years ago
Replying to neteler:
The user took this map
and passed it to R, and then redefined:
and then he wrote b as npp.tif using the QGIS plugin manageR.
Maybe there is a precision issue before, somwhere in the R or QGIS part?
It appears so. The inaccuracy is present in the TIFF file; it isn't being introduced by libtiff, GDAL or GRASS.
I'm not quite sure where it would come from; 1/4 is exactly representable in single-precision floating-point. My first guess would be something calculating 360*(1/1440) rather than 360/1440 (1/1440 isn't exactly representable). This can occur if code is compiled with -funsafe-math-optimizations or similar.
In any case, there's still the issue that GRASS takes a lat/lon region which is actually 360.000000000288 degrees across and treats it as if it's 0.000000000288 degrees across. Wrapping coordinates is one thing; wrapping intervals is something else entirely.
Ultimately, you can't take algorithms which work for Euclidian geometry and make them work for spherical geometry with nothing more than a couple of quick hacks. On a related note, I still haven't got anywhere with the "r.resamp.stats: nulls along western boundary" issue reported by Hamish (this should probably be added as a ticket).
comment:4 by , 15 years ago
Component: | default → Raster |
---|---|
Keywords: | r.in.gdal added |
comment:5 by , 9 years ago
Milestone: | 6.4.0 → 6.4.6 |
---|
comment:6 by , 8 years ago
Replying to glynn:
Replying to neteler:
The user took this map
and passed it to R, and then redefined:
and then he wrote b as npp.tif using the QGIS plugin manageR.
Maybe there is a precision issue before, somwhere in the R or QGIS part?
It appears so. The inaccuracy is present in the TIFF file; it isn't being introduced by libtiff, GDAL or GRASS.
The example file npp.tif is no longer existing. The problem can be reproduced with
g.region n=85.000000000114397 \ s=-58.000000000000007 \ e=180.000000000288 \ w=-180 rows=572 cols=1440 -p
In GRASS 7.2, this gives
projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 85N south: 58S west: 180W east: 179:59:59.999999W <-- wrong EW extent of 0.000001 seconds nsres: 0:15 ewres: 0 <-- invalid because of wrong east rows: 572 cols: 1440 cells: 823680
A subsequent g.region -p
gives
ERROR: Syntax error in cell header
In trunk, setting the region as above gives
projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 85N south: 58S west: 180W east: 180:00:00.000001E <-- ok nsres: 0:15 ewres: 0:15 <-- ok rows: 572 cols: 1440 cells: 823680
The EW extent is now slightly (0.000001 seconds) larger than 360 degrees, but g.region does not complain and east is easy to fix.
Replying to neteler:
It's that extra 0.0000000000002 that's causing the problems:
Note the "east = 180.000000000288" part. The cellhd structure never changes, but the value is wrapped when written to the cellhd file. Essentially, the map is ever so slightly more than 360 degrees wide.
No, it's the wrap-around which causes the problem.
[And which will forever cause problems. We'll just have to work around them as we find them. Almost everything that you "know" about coordinate geometry is wrong for spherical geometry. E.g. just because x1 == x2, that doesn't mean that x1-x0 == x2-x0.]
So how do you want to "fix" this? Clamp east to west+360? Clamp west to east-360? Clamp both to (west+east)/2 +/- 180? What if the map is significantly more than 360 degrees wide?