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

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 :slight_smile: 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.

Hi Marcello,

On Sun, Jul 1, 2012 at 9:17 AM, Marcello Gorini <gorini@gmail.com> wrote:
...

But when trying to do g.region I get

ERROR: Invalid region: Invalid coordinates

Please post

g.proj -p

...

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

...

| Projection: x,y

... this looks suspicious! Above output may help us to find the
problem.

Markus

Markus:

...
> |
Projection: x,y

... this looks suspicious! Above output may help us to find
the problem.

probably he was following my tutorial at
http://grass.osgeo.org/wiki/Global_datasets#ETOPO1_.28DEM.29

which uses a nasty hack to get around the original grid-registered
data exceeding the limits of polar coordinate space:
r.in.bin -f in=etopo1_bed_g.flt out=etopo1_bed_g \
   n=90.008333333335 s=-90.008333333335 e=180.00833333334 \
   w=-180.00833333334 rows=10801 cols=21601 anull=-9999

but the cellhd/ file was not reset in sync with the location.

I have now rewriten that wiki entry to use a r.region cleanup
method instead of ugly filesystem edits.

Hamish

On Sun, Jul 1, 2012 at 9:17 AM, Marcello Gorini <gorini@gmail.com> wrote:

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.

What are the extents of the vector in latlon before projecting into
the mercator location? Note that regions spanning 180 degrees
longitude or more can not be reprojected to mercator locations due to
mathematical constraints.

[...]

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

This inf looks very suspicious, the latlon vector extents very
probably too large and could not be represented in the projection of
the target location.

Markus M

|
| 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 :slight_smile: 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.

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Thanks for the suggestions.

Markus:

Please post
g.proj -p

-PROJ_INFO--------------------

name : Mercator
proj : merc
datum : wgs84
ellps : wgs84
lon_0 : 0
k : 1
x_0 : 0
y_0 : 0
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1

Doesn’t seem to have a problem, does it?

Now I looked at the revised wiki (thanks Hamish) and followed it entirely.

It worked. I guess either I hadn’t followed everything the last time or this r.region part was
really needed.

OK, but I have another problem.

When I first imported the data to the latlog location, besides being flipped, it appeared stretched in the N-S direction.
When I fixed it with region w=-180 e=180, everything fell back to its correct place and it was fine, with no more stretching.

Now, when I reprojected it, it is stretched again.

Doing g.region in my latlong location, I get a resolution of 5 minutes with rows=2159 cols=4321.
When I r.proj -g in the mercator location, following Hamish suggestion, I get
Input Projection Parameters: +proj=longlat +no_defs +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000
Input Unit Factor: 1
Output Projection Parameters: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +no_defs +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000
Output Unit Factor: 1
Input map rbd@ATLANTIC in location <GEBCO_GEODETIC>:
n=46050366.66617828 s=-46050366.66617872 w=-20028233.86541469 e=20028233.86541469 rows=2159 cols=4321

But when I apply the n-s w-e it suggests and r.proj again:

Input:
Cols: 4321 (4321)
Rows: 2159 (2159)
North: 89.916705 (89.916705)
South: -89.916705 (-89.916705)
West: -179.916686 (-179.916686)
East: 179.916686 (179.916686)
EW-res: 0.083275
NS-res: 0.083295

Output:
Cols: 4321 (4321)
Rows: 2159 (2159)
North: 46050366.666178 (46050366.666178)
South: -46050366.666179 (-46050366.666179)
West: -20028233.865415 (-20028233.865415)
East: 20028233.865415 (20028233.865415)
EW-res: 9270.184617
NS-res: 42658.977921

You see? It stretched the NS-res. Is it me or the data? (probably me, of course).

Thanks again,
Marcello.

On Mon, Jul 2, 2012 at 2:47 AM, Hamish <hamish_b@yahoo.com> wrote:

Markus:

|
Projection: x,y

… this looks suspicious! Above output may help us to find
the problem.

probably he was following my tutorial at
http://grass.osgeo.org/wiki/Global_datasets#ETOPO1_.28DEM.29

which uses a nasty hack to get around the original grid-registered
data exceeding the limits of polar coordinate space:
r.in.bin -f in=etopo1_bed_g.flt out=etopo1_bed_g
n=90.008333333335 s=-90.008333333335 e=180.00833333334
w=-180.00833333334 rows=10801 cols=21601 anull=-9999

but the cellhd/ file was not reset in sync with the location.

I have now rewriten that wiki entry to use a r.region cleanup
method instead of ugly filesystem edits.

Hamish

Makus:

What are the extents of the vector in latlon before projecting into
the mercator location? Note that regions spanning 180 degrees
longitude or more can not be reprojected to mercator locations due to
mathematical constraints.20037508.34278924

That’s the result of v.info of the of the result of v.in.region

| N: 89:55:00.138825N S: 89:55:00.138825S |
| E: 179:55:00.069428E W: 179:55:00.069428W

This inf looks very suspicious, the latlon vector extents very
probably too large and could not be represented in the projection of
the target location.

Too much for mercator? I am really lost here…

Thanks.
Marcello.

On Mon, Jul 2, 2012 at 9:03 AM, Marcello Gorini <gorini@gmail.com> wrote:

Makus:

What are the extents of the vector in latlon before projecting into
the mercator location? Note that regions spanning 180 degrees
longitude or more can not be reprojected to mercator locations due to
mathematical constraints.20037508.34278924

That's the result of v.info of the of the result of v.in.region

| N: 89:55:00.138825N S: 89:55:00.138825S
|
| E: 179:55:00.069428E W: 179:55:00.069428W

This inf looks very suspicious, the latlon vector extents very
probably too large and could not be represented in the projection of
the target location.

Too much for mercator? I am really lost here....

I think for mercator this is fine, but not for transverse mercator.

Markus M

Projection: x,y

It is impossible to reproject from a simple XY location to a
projected location, or vice versa. Simple XY is just like graph
paper, with no Earth-based geo-* part to it. It is most commonly
used for imagery where x,y are measured in pixels, and manual geo-
referencing must be performed, not a reprojection between known
coord systems.

Start over from the beginning, creating a new lat/lon location and
importing your data into it. (assuming that's what the source
data's native coordinates are in)

sorry for introducing the nasty hack which suggested using xy,
Hamish

Hamish said:

It is impossible to reproject from a simple XY location to a
projected location, or vice versa. Simple XY is just like graph
paper, with no Earth-based geo-* part to it. It is most commonly
used for imagery where x,y are measured in pixels, and manual geo-
referencing must be performed, not a reprojection between known
coord systems.

Yes, I hadn’t noticed this xy. Sure it won’t work.

Start over from the beginning, creating a new lat/lon location and
importing your data into it. (assuming that’s what the source
data’s native coordinates are in)

I did and it worked flawlessly.

sorry for introducing the nasty hack which suggested using xy,
Hamish

I really don’t remember doing that, but anyway, since you updated
the wiki page, I cannot press charges against you anymore.

Thanks a lot!!

All the best,
Marcello.