Dear all,<br><br>After importing a netCDF file and correcting its bounds using <br>r.region so that it spans from -180/180 instead of 0/360, a tried to<br>reproject it in the way I always do, following the GRASS book.<br><br>
I use v.in.region to get the region as a vector, v.proj in a mercator<br>location and then g.region to the projected vector, before applying <br>r.proj to my raster.<br><br>But when trying to do g.region I get<br><br>ERROR: Invalid region: Invalid coordinates<br>
<br>Following a recent Hamish reply, I noticed that summing the cells<br>I get more than 180 N-S and more than 360 E-W. This time, it's not <br>because the cell size is greater than it should be, but I think that's <br>
because the grid is gridline-registered and than GRASS is "creating"<br>an extra cell in both directions.<br><br>Although this is not a ETOPO1 grid (but it's in a way based on it), I<br>followed GRASS Wiki and did something like this<br>
<br><pre># 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</pre><br>Now, the new grid is within accepted bounds, but the error persists and<br>the result of <a href="http://v.info">v.info</a> of the projected region is<br><br>+----------------------------------------------------------------------------+<br>
| Layer: rbd <br> | Mapset: BRASIL <br> | Location: GEBCO_MERCATOR <br>
| Database: /home/marcello/grassdata <br> | Title: <br> | Map scale: 1:1 <br>
| Map format: native <br> | Name of creator: marcello <br> | Organization: <br>
| Source date: Sun Jul 1 07:55:56 2012 <br> |----------------------------------------------------------------------------<br> | Type of Map: vector (level: 2) <br>
| <br> | Number of points: 0 Number of areas: 0 <br> | Number of lines: 0 Number of islands: 0 <br>
| Number of boundaries: 1 Number of faces: 0 <br> | Number of centroids: 1 Number of kernels: 0 <br> | <br>
| Map is 3D: No <br> | Number of dblinks: 0 <br> | <br>
| Projection: x,y <br> | N: inf S: -4605.12989327 <br> | E: inf W: -20037508.34278924 <br>
| <br> | Digitization threshold: 0 <br> | Comments: <br>
| <br> +----------------------------------------------------------------------------+<br><br>Obviously, something is wrong with the N and E coordinates.<br>
<br>Any ideas?<br><br>Thanks in advance.<br><br>Marcello.<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>the problem is with broken cellsize value, resulting in a<br>
latitude beyond the north pole.<br>
<br>
>> 0.0083333337679505 * 18000 - 60<br>
ans = 90.000007823109<br>
<br>
which is > 90.<br>
<br>
it seems that whatever software exported it (don't be afraid to<br>
name names :) was holding or calculating the resolution with<br>
single-precision floating point numbers but exporting it as if<br>
it were a double-precision number. So the second half the number<br>
is inexact jibberish.<br>
<br>
Edit the cell size back to 0.008333333 (no small feat with a<br>
3.7gb file, even for vi) and it'll work.<br>
<br>
<br>
Importing with GDAL would give the same illegal-north latitude<br>
result, although r.in.gdal now has a '-l' flag to reset the<br>
northern boundary into something legal (after which you Must<br>
repair it to the real value with r.region), although you could<br>
get the same effect by editing the header to lie about the<br>
cellsize or southern value to get it to fit into legal lat/lon,<br>
then again use r.region to set the bounds exactly. (the<br>
resolution as seen with <a href="http://r.info/" target="_blank">r.info</a> should end up exactly at 30 arc-<br>
sec; bypass the built in failsafe checks at your own risk)<br>
But the real solution is to get the software that created it<br>
to not export broken files.<br>