r.in.gdal region troubles in LL

[neteler]

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

[glynn]

Replying to neteler:

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

Pixel Size = (0.250000000000200,-0.250000000000200)

It's that extra 0.0000000000002 that's causing the problems:

> print adfGeoTransform
$2 = {-180, 0.25000000000020001, 0, 85.000000000114397, 0, 
  -0.25000000000020001}

> print cellhd
$3 = {format = -1215084252, compressed = -1208408148, rows = 572, 
  rows3 = 572, cols = 1440, cols3 = 1440, depths = 1, proj = -1215649599, 
  zone = -1215082992, ew_res = 0.25000000000020001, 
  ew_res3 = 0.25000000000020001, ns_res = 0.25000000000020001, 
  ns_res3 = 0.25000000000020001, tb_res = 1, north = 85.000000000114397, 
  south = -58.000000000000007, east = 180.000000000288, west = -180, top = 1, 
  bottom = 0}

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.

Apparently the global wrap-around recognition is failing for 180W/E.

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?

[neteler]

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

[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.

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).

[mmetz]

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.