[GRASS-user] Invalid region and coordinates when trying to reproject.

[Marcello Gorini]

Dear all,

After importing a netCDF file and correcting its bounds using
r.region so that it spans from -180/180 instead of 0/360, a tried to
reproject it in the way I always do, following the GRASS book.

I use v.in.region to get the region as a vector, v.proj in a mercator
location and then g.region to the projected vector, before applying
r.proj to my raster.

But when trying to do g.region I get

ERROR: Invalid region: Invalid coordinates

Following a recent Hamish reply, I noticed that summing the cells
I get more than 180 N-S and more than 360 E-W. This time, it's not
because the cell size is greater than it should be, but I think that's
because the grid is gridline-registered and than GRASS is "creating"
an extra cell in both directions.

Although this is not a ETOPO1 grid (but it's in a way based on it), I
followed GRASS Wiki and did something like this

# reduce region by 1 cell
g.region rast=etopo1_bed_g
eval `g.region -g`
g.region n=n-$nsres s=s+$nsres e=e-$ewres -p

# save smaller raster and remove original
r.mapcalc "etopo1_bed_g.crop = etopo1_bed_g"
g.remove etopo1_bed_g


Now, the new grid is within accepted bounds, but the error persists and
the result of v.info of the projected region is

+----------------------------------------------------------------------------+
 | Layer:           rbd
 | Mapset:
BRASIL
 | Location:
GEBCO_MERCATOR
 | Database:
/home/marcello/grassdata
 |
Title:
 | Map scale:
1:1
 | Map format:      native
 | Name of creator:
marcello
 | Organization:
 | Source date:     Sun Jul  1 07:55:56 2012
 |----------------------------------------------------------------------------
 |   Type of Map:  vector (level: 2)
 |

 |   Number of points:       0               Number of areas:      0
 |   Number of lines:        0               Number of islands:
0
 |   Number of boundaries:   1               Number of faces:      0
 |   Number of centroids:    1               Number of kernels:    0
 |

 |   Map is 3D:              No
 |   Number of dblinks:      0
 |
 |         Projection: x,y
 |               N:               inf    S:
-4605.12989327
 |               E:               inf    W: -20037508.34278924
 |

 |   Digitization threshold:
0
 |   Comments:
 |
 +----------------------------------------------------------------------------+

Obviously, something is wrong with the N and E coordinates.

Any ideas?

Thanks in advance.

Marcello.

















the problem is with broken cellsize value, resulting in a
latitude beyond the north pole.

>> 0.0083333337679505 * 18000 - 60
ans =      90.000007823109

which is > 90.

it seems that whatever software exported it (don't be afraid to
name names :) was holding or calculating the resolution with
single-precision floating point numbers but exporting it as if
it were a double-precision number. So the second half the number
is inexact jibberish.

Edit the cell size back to 0.008333333 (no small feat with a
3.7gb file, even for vi) and it'll work.


Importing with GDAL would give the same illegal-north latitude
result, although r.in.gdal now has a '-l' flag to reset the
northern boundary into something legal (after which you Must
repair it to the real value with r.region), although you could
get the same effect by editing the header to lie about the
cellsize or southern value to get it to fit into legal lat/lon,
then again use r.region to set the bounds exactly. (the
resolution as seen with r.info should end up exactly at 30 arc-
sec; bypass the built in failsafe checks at your own risk)
But the real solution is to get the software that created it
to not export broken files.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20120701/8be4e4d5/attachment-0001.html>