[GRASS-dev] [GRASS GIS] #352: r.in.gdal region troubles in LL

GRASS GRASS developer
,
1

#352: r.in.gdal region troubles in LL
-------------------------+--------------------------------------------------
Reporter: neteler | Owner: grass-dev@lists.osgeo.org
     Type: defect | Status: new
Priority: major | Milestone: 6.4.0
Component: default | Version: unspecified
Keywords: | Platform: Unspecified
      Cpu: Unspecified |
-------------------------+--------------------------------------------------
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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/352&gt;
GRASS GIS <http://grass.osgeo.org>

2

#352: r.in.gdal region troubles in LL
--------------------------+-------------------------------------------------
  Reporter: neteler | Owner: grass-dev@lists.osgeo.org
      Type: defect | Status: new
  Priority: major | Milestone: 6.4.0
Component: default | Version: unspecified
Resolution: | Keywords:
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Comment (by glynn):

Replying to [ticket:352 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?

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/352#comment:1&gt;
GRASS GIS <http://grass.osgeo.org>

3

#352: r.in.gdal region troubles in LL
--------------------------+-------------------------------------------------
  Reporter: neteler | Owner: grass-dev@lists.osgeo.org
      Type: defect | Status: new
  Priority: major | Milestone: 6.4.0
Component: default | Version: unspecified
Resolution: | Keywords:
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Changes (by neteler):

* cc: Agustin.Lobo@ija.csic.es (added)

Comment:

Replying to [comment:1 glynn]:
> Replying to [ticket:352 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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/352#comment:2&gt;
GRASS GIS <http://grass.osgeo.org>

4

#352: r.in.gdal region troubles in LL
--------------------------+-------------------------------------------------
  Reporter: neteler | Owner: grass-dev@lists.osgeo.org
      Type: defect | Status: new
  Priority: major | Milestone: 6.4.0
Component: default | Version: unspecified
Resolution: | Keywords:
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Comment (by glynn):

Replying to [comment:2 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).

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/352#comment:3&gt;
GRASS GIS <http://grass.osgeo.org>

5

#352: r.in.gdal region troubles in LL
--------------------------+-------------------------------------------------
  Reporter: neteler | Owner: grass-dev@lists.osgeo.org
      Type: defect | Status: new
  Priority: major | Milestone: 6.4.0
Component: Raster | Version: unspecified
Resolution: | Keywords: r.in.gdal
  Platform: Unspecified | Cpu: Unspecified
--------------------------+-------------------------------------------------
Changes (by martinl):

  * keywords: => r.in.gdal
  * component: default => Raster

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/352#comment:4&gt;
GRASS GIS <http://grass.osgeo.org>