[GRASS-user] r.in.gdal oddity

Brian_Oney · October 30, 2011, 2:58am

Hello List,

I am importing worldclim .asc files from the ipcc4 simulations (http://ccafs-climate.org/). These are latlong files with a global extent. I have compiled grass 6.4.2svn from source on ubuntu 11.10 and hope that I did everything correctly, according to the wiki.

When importing, here is what I do and what happens:

r.in.gdal -o input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\
hccpr_hadcm3_a2a_2020s_prec_2_5min_asc/prec_10.asc" output="prec_10"

r.info prec_1
  +----------------------------------------------------------------------------+
  | Layer: prec_10@data_2020 Date: Sun Oct 30 02:58:49 2011 |
  | Mapset: data_2020 Login of Creator: root |
  | Location: worldclim |
  | DataBase: /D/GIS_Dbase |
  | Title: ( prec_10 ) |
  | Timestamp: none |
  |----------------------------------------------------------------------------|
  | |
  | Type of Map: raster Number of Categories: 1185 |
  | Data Type: CELL |
  | Rows: 3600 |
  | Columns: 8640 |
  | Total Cells: 31104000 |
  | Projection: Latitude-Longitude |
  | N: 90N S: 60S Res: 0:02:30 |
  | E: 179:59:59.999999W W: 180W Res: 0 |
  | Range of data: min = 0 max = 1185 |
  | |
  | Data Description: |
  | generated by r.in.gdal |
  | |
  | Comments: |
  | r.in.gdal -o input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\ |
  | hccpr_hadcm3_a2a_2020s_prec_2_5min_asc/prec_10.asc" output="prec_10" |
  | |
  +----------------------------------------------------------------------------+

### The resolution is obviously the problem here and I am not sure why it does this. When I try to ### force it into the region with the -l flag I get the same result.
### Then I tell it how it should appear. This corresponds to the regions i.e. g.region -g
r.region map=prec_10 n=90 s=-60 e=180 w=-180
### Have a look
r.info prec_10@data_2020

WARNING: Unable to read header file for raster map
<prec_10@data_2020@data_2020>. Invalid format.line 9: <e-w resol:
          0>
  +----------------------------------------------------------------------------+
  | Layer: prec_10@data_2020 Date: Sun Oct 30 01:14:45 2011 |
  | Mapset: data_2020 Login of Creator: root |
  | Location: worldclim |
  | DataBase: /D/GIS_Dbase |
  | Title: ( prec_10 ) |
  | Timestamp: none |
  |----------------------------------------------------------------------------|
  | |
  | Type of Map: raster Number of Categories: 1185 |
  | Data Type: CELL |
  | |
  | Data Description: |
  | generated by r.in.gdal |
  | |
  | Comments: |
  | r.in.gdal -o -l input="/D/Documents/Spatial_Datasets/Climate/WorldCl\ |
  | im/hccpr_hadcm3_a2a_2020s_prec_2_5min_asc/prec_10.asc" output="prec_\ |
  | 10" |
  | |
  +----------------------------------------------------------------------------+

I have tried the align option with r.region and I am stuck. I have already perused the documentation concerning the .bil file import... I could try that but this should work too. The r.in.arc module gives me similar results.
I really appreciate the help.

Brian

hellik · October 30, 2011, 12:29pm

[...]

I have tried the align option with r.region and I am stuck. I have
already perused the documentation concerning the .bil file import... I
could try that but this should work too. The r.in.arc module gives me
similar results.
I really appreciate the help.

for further information see:

http://osgeo-org.1803224.n2.nabble.com/r-in-gdal-worldclim-data-G-set-window-geographic-latitude-not-valid-for-north-td5975048.html

or

http://pvanb.wordpress.com/2010/11/30/importing-worldclim-climate-bil-datalayers-in-grass-gis-ii/

Helmut

--
View this message in context: http://osgeo-org.1803224.n2.nabble.com/r-in-gdal-oddity-tp6944754p6945447.html
Sent from the Grass - Users mailing list archive at Nabble.com.

hellik · October 30, 2011, 1:41pm

I am importing worldclim .asc files from the ipcc4 simulations
(http://ccafs-climate.org/). These are latlong files with a global
extent. I have compiled grass 6.4.2svn from source on ubuntu 11.10 and
hope that I did everything correctly, according to the wiki.

When importing, here is what I do and what happens:

r.in.gdal -o input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\
hccpr_hadcm3_a2a_2020s_prec_2_5min_asc/prec_10.asc" output="prec_10"

tested here with
http://ccafs-climate.org/data/A2a_2020s/hccpr_hadcm3/hccpr_hadcm3_a2a_2020s_prec_2_5min_asc.zip

gdalinfo hccpr_hadcm3_a2a_2020s_prec_2_5min.asc

Driver: AAIGrid/Arc/Info ASCII Grid
Files: prec_10.asc
Size is 8640, 3600
Coordinate System is `'
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000)
Lower Left (-180.0000000, -60.0000000)
Upper Right ( 180.0000000, 90.0000000)
Lower Right ( 180.0000000, -60.0000000)
Center ( 0.0000000, 15.0000000)
Band 1 Block=8640x1 Type=Int32, ColorInterp=Undefined
NoData Value=-9999

Helmut

--
View this message in context: http://osgeo-org.1803224.n2.nabble.com/r-in-gdal-oddity-tp6944754p6945538.html
Sent from the Grass - Users mailing list archive at Nabble.com.

Brian_Oney · October 30, 2011, 10:25pm

Hi Helmut,

I have read that and done that. I get the same information. Yes, in this case I am supposed to use the -l flag, but I get the same result (as I wrote).

A bit more peculiar when importing the ArcGRID files...

r.info prec_1
  +----------------------------------------------------------------------------+
  | Layer: prec_1 Date: Sun Oct 30 22:57:26 2011 |
  | Mapset: data_1975 Login of Creator: root |
  | Location: worldclim |
  | DataBase: /D/GIS_Dbase |
  | Title: ( prec_1 ) |
  | Timestamp: none |
  |----------------------------------------------------------------------------|
  | |
  | Type of Map: raster Number of Categories: 916 |
  | Data Type: CELL |
  | Rows: 3600 |
  | Columns: 8640 |
  | Total Cells: 31104000 |
  | Projection: Latitude-Longitude |
  | N: 90N S: 60S Res: 0:02:30 |
  | E: 179:59:59.932408W W: 180W Res: 0:00:00.000008 |
  | Range of data: min = 0 max = 916 |
  | |
  | Data Description: |
  | generated by r.in.gdal |
  | |
  | Comments: |
  | r.in.gdal -l input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\ |
  | prec_2-5m_esri/prec/prec_1" output="prec_1" |
  | |
  +----------------------------------------------------------------------------+
## Same thing
r.region map=prec_1 n=90 s=-60 e=180 w=-180

r.info prec_1
  +----------------------------------------------------------------------------+
  | Layer: prec_1 Date: Sun Oct 30 22:57:26 2011 |
  | Mapset: data_1975 Login of Creator: root |
  | Location: worldclim |
  | DataBase: /D/GIS_Dbase |
  | Title: ( prec_1 ) |
  | Timestamp: none |
  |----------------------------------------------------------------------------|
  | |
  | Type of Map: raster Number of Categories: 916 |
  | Data Type: CELL |
  | Rows: 3600 |
  | Columns: 8640 |
  | Total Cells: 31104000 |
  | Projection: Latitude-Longitude |
  | N: 90N S: 60S Res: 0:02:30 |
  | E: 180E W: 180W Res: 0:02:30 |
  | Range of data: min = 0 max = 916 |
  | |
  | Data Description: |
  | generated by r.in.gdal |
  | |
  | Comments: |
  | r.in.gdal -l input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\ |
  | prec_2-5m_esri/prec/prec_1" output="prec_1" |
  | |
  +----------------------------------------------------------------------------+
Now r.region lets me pull and stretch as I wish. Notice the slight difference in e-w resolution. Maybe that could be the key. Nonetheless, it is odd that even with overriding the projection check AND OR forcing the map to fit within the bounds, the map is still a bit "skinny"...
I realize there is an Illegal latitude. I have been playing with this for almost a week now... I have RTFM plenty of times. Not to say that I didn't or couldn't miss something, but it would less likely after such a long time of tinkering :). The more pragmatic solution would be to download the .bil files, because I have already worked with them... I can get the thing to work another way, however I want to understand what is wrong and where.

Still stumped. A bit of context may help: I am a semi-experienced grass user (+1.5 years) but just switched to linux. Any common experiences when starting anew? Could I have compiled something incorrectly? It is a bit odd. Maybe you could try importing and tell me what it says? r.in.arc gives me the same result...

Thanks for the help!

Cheers,
Brian

On 10/30/2011 02:41 PM, Helmut Kudrnovsky wrote:

I am importing worldclim .asc files from the ipcc4 simulations
(http://ccafs-climate.org/). These are latlong files with a global
extent. I have compiled grass 6.4.2svn from source on ubuntu 11.10 and
hope that I did everything correctly, according to the wiki.

When importing, here is what I do and what happens:

r.in.gdal -o input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\
hccpr_hadcm3_a2a_2020s_prec_2_5min_asc/prec_10.asc" output="prec_10"

tested here with
http://ccafs-climate.org/data/A2a_2020s/hccpr_hadcm3/hccpr_hadcm3_a2a_2020s_prec_2_5min_asc.zip

gdalinfo hccpr_hadcm3_a2a_2020s_prec_2_5min.asc

Driver: AAIGrid/Arc/Info ASCII Grid
Files: prec_10.asc
Size is 8640, 3600
Coordinate System is `'
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000)
Lower Left (-180.0000000, -60.0000000)
Upper Right ( 180.0000000, 90.0000000)
Lower Right ( 180.0000000, -60.0000000)
Center ( 0.0000000, 15.0000000)
Band 1 Block=8640x1 Type=Int32, ColorInterp=Undefined
NoData Value=-9999

Helmut

--
View this message in context: http://osgeo-org.1803224.n2.nabble.com/r-in-gdal-oddity-tp6944754p6945538.html
Sent from the Grass - Users mailing list archive at Nabble.com.
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

hellik · October 30, 2011, 11:27pm

I have read that and done that. I get the same information. Yes, in this
case I am supposed to use the -l flag, but I get the same result (as I
wrote).
[...]
Now r.region lets me pull and stretch as I wish. Notice the slight
difference in e-w resolution. Maybe that could be the key.

AFAIK there seems to be some precision issue (e-w resolution)
in the producing work flow of this data. I've opened the asc and the ESRI
grid also in
ArcGIS 10, there is this same little difference in the data extent.

Nonetheless,
it is odd that even with overriding the projection check AND OR forcing
the map to fit within the bounds, the map is still a bit "skinny"...

you're right. maybe Markus M (cc'ed) has some more insight into this?

Helmut

--
View this message in context: http://osgeo-org.1803224.n2.nabble.com/r-in-gdal-oddity-tp6944754p6946588.html
Sent from the Grass - Users mailing list archive at Nabble.com.

Markus_GRASS · October 31, 2011, 8:15am

Helmut Kudrnovsky wrote:

I am importing worldclim .asc files from the ipcc4 simulations
(http://ccafs-climate.org/). These are latlong files with a global
extent. I have compiled grass 6.4.2svn from source on ubuntu 11.10 and
hope that I did everything correctly, according to the wiki.

When importing, here is what I do and what happens:

r.in.gdal -o input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\
hccpr_hadcm3_a2a_2020s_prec_2_5min_asc/prec_10.asc" output="prec_10"

tested here with
http://ccafs-climate.org/data/A2a_2020s/hccpr_hadcm3/hccpr_hadcm3_a2a_2020s_prec_2_5min_asc.zip

gdalinfo hccpr_hadcm3_a2a_2020s_prec_2_5min.asc

Driver: AAIGrid/Arc/Info ASCII Grid
Files: prec_10.asc
Size is 8640, 3600
Coordinate System is `'
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000)
Lower Left (-180.0000000, -60.0000000)
Upper Right ( 180.0000000, 90.0000000)
Lower Right ( 180.0000000, -60.0000000)
Center ( 0.0000000, 15.0000000)
Band 1 Block=8640x1 Type=Int32, ColorInterp=Undefined
NoData Value=-9999

There is a tiny precision loss:
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
should be
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.041666666666667,-0.041666666666667)

pixel size is now as close as possible to 2.5 min.

This is neither a gdal nor a grass bug, but a bug of the software used
to produce these data.

This can easily be fixed with e.g.
gdal_translate -a_ullr -180 90 180 -60 prec_1.asc prec_1.tif

then import the corrected tif.

HTH,

Markus M

Brian_Oney · October 31, 2011, 7:42pm

Hi Markus M,
Thanks, it works. Cooking with gas now.

It is still odd though. It is not just an imprecision I am afraid. When I said "skinny", I was being modest.
If you look a bit closer, it is not a minor error in import. The map (and the rest as well) is w=-180 e=-179.999... and north and south are normal. It should have a global extent.

A bit of additional info: when the location lacks a projection, I don't have problems.

Anyways, I appreciate the help.

Cheers,
Brian

On 10/31/2011 09:15 AM, Markus Metz wrote:

Helmut Kudrnovsky wrote:

I am importing worldclim .asc files from the ipcc4 simulations
(http://ccafs-climate.org/). These are latlong files with a global
extent. I have compiled grass 6.4.2svn from source on ubuntu 11.10 and
hope that I did everything correctly, according to the wiki.

When importing, here is what I do and what happens:

r.in.gdal -o input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\
hccpr_hadcm3_a2a_2020s_prec_2_5min_asc/prec_10.asc" output="prec_10"

tested here with
http://ccafs-climate.org/data/A2a_2020s/hccpr_hadcm3/hccpr_hadcm3_a2a_2020s_prec_2_5min_asc.zip

gdalinfo hccpr_hadcm3_a2a_2020s_prec_2_5min.asc

Driver: AAIGrid/Arc/Info ASCII Grid
Files: prec_10.asc
Size is 8640, 3600
Coordinate System is `'
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000)
Lower Left (-180.0000000, -60.0000000)
Upper Right ( 180.0000000, 90.0000000)
Lower Right ( 180.0000000, -60.0000000)
Center ( 0.0000000, 15.0000000)
Band 1 Block=8640x1 Type=Int32, ColorInterp=Undefined
NoData Value=-9999

There is a tiny precision loss:
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
should be
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.041666666666667,-0.041666666666667)

pixel size is now as close as possible to 2.5 min.

This is neither a gdal nor a grass bug, but a bug of the software used
to produce these data.

This can easily be fixed with e.g.
gdal_translate -a_ullr -180 90 180 -60 prec_1.asc prec_1.tif

then import the corrected tif.

HTH,

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

Markus_GRASS · October 31, 2011, 7:59pm

Brian Oney wrote:

Hi Markus M,
Thanks, it works. Cooking with gas now.

So the import with the corrected raster is ok, i.e. is of global (-180
to 180 E, -60 to 90 N) extent?

It is still odd though. It is not just an imprecision I am afraid. When I
said "skinny", I was being modest.
If you look a bit closer, it is not a minor error in import. The map (and
the rest as well) is w=-180 e=-179.999... and north and south are normal. It
should have a global extent.

The skinny appearance (just a tiny slice of data) happens because the
input data exceed 360 degree in E-W extent by a tiny fraction. This is
automatically adjusted such that the input data do not exceed 360
degree in E-W direction. Maybe this automated adjustment needs some
improvement?

A bit of additional info: when the location lacks a projection, I don't have
problems.

As expected, it's a latlon issue.

Markus M

Anyways, I appreciate the help.

Cheers,
Brian

On 10/31/2011 09:15 AM, Markus Metz wrote:

Helmut Kudrnovsky wrote:

I am importing worldclim .asc files from the ipcc4 simulations
(http://ccafs-climate.org/). These are latlong files with a global
extent. I have compiled grass 6.4.2svn from source on ubuntu 11.10 and
hope that I did everything correctly, according to the wiki.

When importing, here is what I do and what happens:

r.in.gdal -o input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\
hccpr_hadcm3_a2a_2020s_prec_2_5min_asc/prec_10.asc" output="prec_10"

tested here with

http://ccafs-climate.org/data/A2a_2020s/hccpr_hadcm3/hccpr_hadcm3_a2a_2020s_prec_2_5min_asc.zip

gdalinfo hccpr_hadcm3_a2a_2020s_prec_2_5min.asc

Driver: AAIGrid/Arc/Info ASCII Grid
Files: prec_10.asc
Size is 8640, 3600
Coordinate System is `'
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000)
Lower Left (-180.0000000, -60.0000000)
Upper Right ( 180.0000000, 90.0000000)
Lower Right ( 180.0000000, -60.0000000)
Center ( 0.0000000, 15.0000000)
Band 1 Block=8640x1 Type=Int32, ColorInterp=Undefined
NoData Value=-9999

There is a tiny precision loss:
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
should be
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.041666666666667,-0.041666666666667)

pixel size is now as close as possible to 2.5 min.

This is neither a gdal nor a grass bug, but a bug of the software used
to produce these data.

This can easily be fixed with e.g.
gdal_translate -a_ullr -180 90 180 -60 prec_1.asc prec_1.tif

then import the corrected tif.

HTH,

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

Brian_Oney · October 31, 2011, 8:25pm

Yes, everything works great. Thanks again.

Ok, good to know. Yeah, the algorithm could be improved I guess. But, the improvement that I would suggest is supposed to be the -l flag, right? Just force it to be 180W and 180E... I thought (the -l flag) it is just a recent addition to GRASS 6.4, so maybe it needs a some t.l.c...

Thanks again.

Cheers,
Brian

On 10/31/2011 08:59 PM, Markus Metz wrote:

Brian Oney wrote:

Hi Markus M,
Thanks, it works. Cooking with gas now.

So the import with the corrected raster is ok, i.e. is of global (-180
to 180 E, -60 to 90 N) extent?

It is still odd though. It is not just an imprecision I am afraid. When I
said "skinny", I was being modest.
If you look a bit closer, it is not a minor error in import. The map (and
the rest as well) is w=-180 e=-179.999... and north and south are normal. It
should have a global extent.

The skinny appearance (just a tiny slice of data) happens because the
input data exceed 360 degree in E-W extent by a tiny fraction. This is
automatically adjusted such that the input data do not exceed 360
degree in E-W direction. Maybe this automated adjustment needs some
improvement?

A bit of additional info: when the location lacks a projection, I don't have
problems.

As expected, it's a latlon issue.

Markus M

Anyways, I appreciate the help.

Cheers,
Brian

On 10/31/2011 09:15 AM, Markus Metz wrote:

Helmut Kudrnovsky wrote:

I am importing worldclim .asc files from the ipcc4 simulations
(http://ccafs-climate.org/). These are latlong files with a global
extent. I have compiled grass 6.4.2svn from source on ubuntu 11.10 and
hope that I did everything correctly, according to the wiki.

When importing, here is what I do and what happens:

r.in.gdal -o input="/D/Documents/Spatial_Datasets/Climate/WorldClim/\
hccpr_hadcm3_a2a_2020s_prec_2_5min_asc/prec_10.asc" output="prec_10"

tested here with

http://ccafs-climate.org/data/A2a_2020s/hccpr_hadcm3/hccpr_hadcm3_a2a_2020s_prec_2_5min_asc.zip

gdalinfo hccpr_hadcm3_a2a_2020s_prec_2_5min.asc

Driver: AAIGrid/Arc/Info ASCII Grid
Files: prec_10.asc
Size is 8640, 3600
Coordinate System is `'
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000)
Lower Left (-180.0000000, -60.0000000)
Upper Right ( 180.0000000, 90.0000000)
Lower Right ( 180.0000000, -60.0000000)
Center ( 0.0000000, 15.0000000)
Band 1 Block=8640x1 Type=Int32, ColorInterp=Undefined
NoData Value=-9999

There is a tiny precision loss:
Origin = (-180.000000000000000,90.000000000119996)
Pixel Size = (0.041666666666700,-0.041666666666700)
should be
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.041666666666667,-0.041666666666667)

pixel size is now as close as possible to 2.5 min.

This is neither a gdal nor a grass bug, but a bug of the software used
to produce these data.

This can easily be fixed with e.g.
gdal_translate -a_ullr -180 90 180 -60 prec_1.asc prec_1.tif

then import the corrected tif.

HTH,

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

neteler · November 1, 2011, 1:46pm

On Mon, Oct 31, 2011 at 9:25 PM, Brian Oney <zenlines@gmail.com> wrote:

Yes, everything works great. Thanks again.

Ok, good to know. Yeah, the algorithm could be improved I guess. But, the
improvement that I would suggest is supposed to be the -l flag, right?

Yes, that was the idea of the -l flag.

Just force it to be 180W and 180E... I thought (the -l flag) it is just a recent
addition to GRASS 6.4, so maybe it needs a some t.l.c...

I run into the same problem recently (see my other posting) and didn't
manage to import with -l. Using

gdal_translate -a_ullr -180 90 180 -60 inmap outmap

is a solution but it would be better to get -l working (since it should do
the same). The code piece is

    /* constrain to geographic coords */
    if (flag_l->answer && G_projection() == PROJECTION_LL) {
        if (cellhd.north > 90.) cellhd.north = 90.;
        if (cellhd.south < -90.) cellhd.south = -90.;
        if (cellhd.east > 360.) cellhd.east = 180.;
        if (cellhd.west < -180.) cellhd.west = -180.;
        cellhd.ns_res = (cellhd.north - cellhd.south) / cellhd.rows;
        cellhd.ew_res = (cellhd.east - cellhd.west) / cellhd.cols;
        cellhd.ew_res3 = cellhd.ew_res;
        cellhd.ns_res3 = cellhd.ns_res;

        G_warning(_("Map bounds have been constrained to geographic "
            "coordinates. You will almost certainly want to check "
            "map bounds and resolution with r.info and reset them "
            "with r.region before going any further."));
    }

... the warning is not particularly clear to me - "reset them" in this sense?

Markus

Markus_GRASS · November 1, 2011, 2:33pm

On Tue, Nov 1, 2011 at 2:46 PM, Markus Neteler <neteler@osgeo.org> wrote:

On Mon, Oct 31, 2011 at 9:25 PM, Brian Oney <zenlines@gmail.com> wrote:

Yes, everything works great. Thanks again.

Ok, good to know. Yeah, the algorithm could be improved I guess. But, the
improvement that I would suggest is supposed to be the -l flag, right?

Yes, that was the idea of the -l flag.

Just force it to be 180W and 180E... I thought (the -l flag) it is just a recent
addition to GRASS 6.4, so maybe it needs a some t.l.c...

I run into the same problem recently (see my other posting) and didn't
manage to import with -l. Using

gdal_translate -a_ullr -180 90 180 -60 inmap outmap

is a solution but it would be better to get -l working (since it should do
the same). The code piece is

/* constrain to geographic coords */
if (flag_l->answer && G_projection() == PROJECTION_LL) {
if (cellhd.north > 90.) cellhd.north = 90.;
if (cellhd.south < -90.) cellhd.south = -90.;

this:

   if \(cellhd\.east &gt; 360\.\) cellhd\.east = 180\.;

should probably be
if (cellhd.east > 180.) cellhd.east = 180.;

But I faintly remember some maps to cover 0, 360 instead of -180, 180.

Furthermore, a map spanning 170,190 (the datumline) should stay as is,
the correction is probably only necessary if E - W > 360, E - W as
provided by GDAL, before GRASS does any conversion.

Also, it may be safer to do the correction after the projection check
and optional creation of a new location, not before.

Markus M

   if \(cellhd\.west &lt; \-180\.\) cellhd\.west = \-180\.;
   cellhd\.ns\_res = \(cellhd\.north \- cellhd\.south\) / cellhd\.rows;
   cellhd\.ew\_res = \(cellhd\.east \- cellhd\.west\) / cellhd\.cols;
   cellhd\.ew\_res3 = cellhd\.ew\_res;
   cellhd\.ns\_res3 = cellhd\.ns\_res;

   G\_warning\(\_\(&quot;Map bounds have been constrained to geographic &quot;
       &quot;coordinates\. You will almost certainly want to check &quot;
       &quot;map bounds and resolution with r\.info and reset them &quot;
       &quot;with r\.region before going any further\.&quot;\)\);

}

... the warning is not particularly clear to me - "reset them" in this sense?

Markus

Brian_Oney · November 1, 2011, 2:45pm

One would think the -l flag would do what it is supposed to, but this then squeezes the raster into a sliver. I get the warning...

May I ask why?:
...
if (cellhd.east > 360.) cellhd.east = 180.;
...
Why not?
...
if (cellhd.east > 180.) cellhd.east = 180.;
...

Also, is the 3d raster stuff necessary?

A further problem appears, when trying to "reset" the raster, which has an e-w resolution of 0 to have a global extent with a corresponding resolution. Maybe the r.region needs a funny something like:

if (e-w_resol == 0.) e-w_resol= 0.000000001; /* or something like that */

Thanks for all the info!

Cheers,
Brian

On 11/01/2011 02:46 PM, Markus Neteler wrote:

On Mon, Oct 31, 2011 at 9:25 PM, Brian Oney<zenlines@gmail.com> wrote:

Yes, everything works great. Thanks again.

Ok, good to know. Yeah, the algorithm could be improved I guess. But, the
improvement that I would suggest is supposed to be the -l flag, right?

Yes, that was the idea of the -l flag.

Just force it to be 180W and 180E... I thought (the -l flag) it is just a recent
addition to GRASS 6.4, so maybe it needs a some t.l.c...

I run into the same problem recently (see my other posting) and didn't
manage to import with -l. Using

  gdal_translate -a_ullr -180 90 180 -60 inmap outmap

is a solution but it would be better to get -l working (since it should do
the same). The code piece is

     /* constrain to geographic coords */
     if (flag_l->answer&& G_projection() == PROJECTION_LL) {
         if (cellhd.north> 90.) cellhd.north = 90.;
         if (cellhd.south< -90.) cellhd.south = -90.;
         if (cellhd.east> 360.) cellhd.east = 180.;
         if (cellhd.west< -180.) cellhd.west = -180.;
         cellhd.ns_res = (cellhd.north - cellhd.south) / cellhd.rows;
         cellhd.ew_res = (cellhd.east - cellhd.west) / cellhd.cols;
         cellhd.ew_res3 = cellhd.ew_res;
         cellhd.ns_res3 = cellhd.ns_res;

         G_warning(_("Map bounds have been constrained to geographic "
             "coordinates. You will almost certainly want to check "
             "map bounds and resolution with r.info and reset them "
             "with r.region before going any further."));
     }

... the warning is not particularly clear to me - "reset them" in this sense?

Markus

Hamish · November 3, 2011, 4:17am

Hi, coming to this thread a bit late, sorry if I miss the point;
I'll study the rest of the thread/situation better asap, but some
general comments--

Markus wrote:

> Just force it to be 180W and 180E... I thought (the -l
> flag) it is just a recent addition to GRASS 6.4, so maybe
> it needs a some t.l.c...

AFAIK it is working as intended and without bugs.

The -l flag is intended to get the raster data into the lat/lon
location, at the expense of breaking (to a greater or lesser
extent) the bounds. When you use that flag, you pretty much Must
run r.region (and/or g.region + r.mapcalc in a crop operation)
before doing anything else. Converting from grid-confluence to
cell-center convention could also require working in a 1/2 cell
shift.

see the "Error Messages" section of the r.in.gdal help page
for expanded explanation and help with the -l flag.

As MarkusM noted, 0-360 is valid and used by planetary scientists,
and I notice some modern GRIB raster data from the US National
Weather Service/US Navy ocean models use that now too.

The -l flag's longitude fix is really just a while-we're-here
thing, the main adjustment the flag is intended to deal with is
the hard +/-90 deg latitude limit.

    G\_warning\(\_\(&quot;Map bounds have been constrained to geographic &quot;
        &quot;coordinates\. You will almost certainly want to check &quot;
        &quot;map bounds and resolution with r\.info and reset them &quot;
        &quot;with r\.region before going any further\.&quot;\)\);
\}

... the warning is not particularly clear to me - "reset
them" in this sense?

"You will almost certainly want to check map bounds and
resolution with r.info and reset them with r.region before
going any further."

maybe I'm not understanding the question or missing an implied
noun, but "them" == map bounds and resolution in this context.

best,
Hamish

neteler · November 3, 2011, 6:21am

On Thu, Nov 3, 2011 at 5:17 AM, Hamish <hamish_b@yahoo.com> wrote:

Hi, coming to this thread a bit late, sorry if I miss the point;
I'll study the rest of the thread/situation better asap, but some
general comments--

Markus wrote:

> Just force it to be 180W and 180E... I thought (the -l
> flag) it is just a recent addition to GRASS 6.4, so maybe
> it needs a some t.l.c...

AFAIK it is working as intended and without bugs.

Please check then
http://article.gmane.org/gmane.comp.gis.grass.user/41407

Here I don't see that -l works.

thanks
Markus

Brian_Oney · November 3, 2011, 10:28am

Hmmm, maybe just a suggestion, although it won't fix the current problem but to make the code better (I hope)... because it apperas the if one has an imprecise 0-360 map that thing could go awry, which may justify another flag (if I may) or something similar to the few changes below.

     /* constrain to geographic coords */
     if (flag_l->answer&& G_projection() == PROJECTION_LL) {
         if (cellhd.north> 90.) cellhd.north = 90.;
         if (cellhd.south< -90.) cellhd.south = -90.;
         if (cellhd.east> 360.) cellhd.east = 360.;

/* for the case of imprecise coordinates 0-360: if the difference is greater than 360 AND it is the "other" ll projection AND the cellhd.west is strange THEN set cellhd.west to 0*/
if (cellhd.east - cellhd.west> 360.&& cellhd.east>= 181.&& cellhd.west< 0.) cellhd.west = 0.;

  /* maybe just one degree of imprecision allowed? */
  if (cellhd.east> 180.&& cellhd.east< 181.) cellhd.east = 180.;
         if (cellhd.west< -180.) cellhd.west = -180.;
         cellhd.ns_res = (cellhd.north - cellhd.south) / cellhd.rows;
         cellhd.ew_res = (cellhd.east - cellhd.west) / cellhd.cols;
         cellhd.ew_res3 = cellhd.ew_res;
         cellhd.ns_res3 = cellhd.ns_res;

         G_warning(_("Map bounds have been constrained to geographic "
             "coordinates. You will almost certainly want to check "
             "map bounds and resolution with r.info and reset them "
             "with r.region before going any further."));
     }

I am not sure how to perform math operations in if-expressions in C, but above is the general idea.

Again, I doubt this will fix the r.in.gdal oddity, but it may improve the -l flag.

Notice that the resolution in the mentioned post is the same as my case... prob just a coincidence.

cheers,
Brian

On 11/03/2011 07:21 AM, Markus Neteler wrote:

On Thu, Nov 3, 2011 at 5:17 AM, Hamish<hamish_b@yahoo.com> wrote:

Hi, coming to this thread a bit late, sorry if I miss the point;
I'll study the rest of the thread/situation better asap, but some
general comments--

Markus wrote:

Just force it to be 180W and 180E... I thought (the -l
flag) it is just a recent addition to GRASS 6.4, so maybe
it needs a some t.l.c...

AFAIK it is working as intended and without bugs.

Please check then
http://article.gmane.org/gmane.comp.gis.grass.user/41407

Here I don't see that -l works.

thanks
Markus

Hamish · November 4, 2011, 6:01am

Hamish:

>> AFAIK it is working as intended and without bugs.

Markus:

> Please check then
> http://article.gmane.org/gmane.comp.gis.grass.user/41407
>
> Here I don't see that -l works.

I will look into it, but my suspicion is that the dataset in question has
IEEE FP issues like -180.00000000000001 which gdalinfo's printf is hiding.

Also to confirm that there hasn't been a new bug introduced I will retry
to import e.g. the ETOPO5 dataset which exists 180W to 180E, and used to
import without any trouble. (I've tested this a lot in the past because I
cross the dateline a lot in my work, and for the last many years it has
worked fine for rasters)

WRT the -l flag, if you use 'r.in.gdal -l' you _must_ use r.region after
to correct the map bounds+resolution. By definition the flag just makes
something up to fool the system into letting you import it, and it is
no surprise that the resulting bounds+resolution are not correct.
maybe we should make rasters imported with -l scale from 0.0 - 1.0 EW,NS
to make it more obvious that the bounds+res are only placeholders, but
I'd guess that would cause more trouble that it fixes.

either way, r.region is the quick-fix solution.

regards,
Hamish

Hamish · November 4, 2011, 6:53am

Markus:

> > Please check then
> > http://article.gmane.org/gmane.comp.gis.grass.user/41407

Hamish:

I will look into it, but my suspicion is that the dataset
in question has IEEE FP issues like -180.00000000000001 which
gdalinfo's printf is hiding.

yup.

the data file's header is broken causing a cumulative rounding problem.
here's the header from the ascii version of it:

ncols 360
nrows 143
xllcorner -180
yllcorner -58
cellsize 1.0000000000008
NODATA_value -9999

change cellsize to 1.000000000000 and then r.in.arc works as expected.
I'd expect the issue to be the same with the binary formats, just not
as easy to repair.

Ben:

It appears that there is a completely way out value in the millions
range, although most of the data is less than a thousand.

it would seem that in places like New York, Tokyo, Mexico City, and
Bangladesh, etc, there are some millions of people who live within the
same rather crowded 100km x 100km area.

# r.univar works on the current computational region (window):
GRASS> g.region rast=gl_gpwfe_pcount_10
GRASS> r.univar -e gl_gpwfe_pcount_10 percentile=99.95
100%
total null and non-null cells: 51480
total null cells: 0

Of the non-null cells:
----------------------
n: 51480
minimum: 0
maximum: 3.16715e+07
range: 3.16715e+07
mean: 132585
mean of absolute values: 132585
standard deviation: 783960
variance: 6.14593e+11
variation coefficient: 591.287 %
sum: 6825499706.1993026733
1st quartile: 0
median (even number of cells): 0
3rd quartile: 1129.64
99.95 percentile: 1.34546e+07

Hamish