[GRASSLIST:4826] Problems importing Data into lat/lon with negtiv south value

GRASS GRASS user
,
1

Dear Users,

I am trying to import some MODIS data into GRASS. In a first step I
managed to convert a subdataset of the MODIS data into img-format using
gdal_translate.

The dataset is a global dataset which should be in lat/long.

gdal_translate HDF4_SDS:EOS_GRID:"modis.hdf":2 -of HFA test_hdf.img
Input file size is 3600, 1800
0...10...20...30...40...50...60...70...80...90...100 - done.
steph@bim:~/modis > gdalinfo test_hdf.img
Driver: HFA/Erdas Imagine Images (.img)
Size is 3600, 1800
Coordinate System is `'
Origin = (0.000000,0.000000)
Pixel Size = (1.00000000,1.00000000)
Metadata:
  HDFEOSVersion=HDFEOS_V2.4
  LOCALGRANULEID=MOD08D3H.A2003090.004.2003093014745.hdf
  PRODUCTIONDATETIME=2003-04-03T01:47:45.000Z
  DAYNIGHTFLAG=Both
  REPROCESSINGACTUAL=processed once
  LOCALVERSIONID=004
  REPROCESSINGPLANNED=further update is anticipated
  VERSIONID=4
  SHORTNAME=MOD08D3H
  INPUTPOINTER=MOD04_L2.A2003090.0250.004.2003092034313.hdf,

[more metadata skipped]

EASTBOUNDINGCOORDINATE=180.000000
WESTBOUNDINGCOORDINATE=-180.000000 SOUTHBOUNDINGCOORDINATE=-90.000000
  NORTHBOUNDINGCOORDINATE=90.000000
  RANGEENDINGDATE=2003-04-01
  RANGEENDINGTIME=00:00:00
  RANGEBEGINNINGDATE=2003-03-31
  RANGEBEGINNINGTIME=00:00:00
  PGEVERSION=4.0.3
  ASSOCIATEDSENSORSHORTNAME=MODIS
  ASSOCIATEDPLATFORMSHORTNAME=Terra
  ASSOCIATEDINSTRUMENTSHORTNAME=MODIS
  PROCESSINGENVIRONMENT=Linux minion425 2.4.18-27.7.xsmp #1 SMP Fri Mar
14 05:52:30 EST 2003 i686 unknown
LOCALINPUTGRANULEID=MOD04_L2.A2003090.0250.004.2003092034313.hdf,

[more metadata skipped]

  LONGNAME=MODIS/Terra Aerosol Cloud Water Vapor Ozone Daily L3 Global
0.1 Deg CMG valid_range=0, 18000
  _FillValue=-9999
  long_name=Sensor Zenith Angle (Cell to Sensor): Mean
  units=Degrees
  scale_factor=0.01
  add_offset=0
  Level_2_Pixel_Values_Read_As=Real
  Derived_From_Level_2_Data_Set=Sensor_Zenith
  Included_Level_2_Nighttime_Data=False
  Quality_Assurance_Data_Set=None
  Statistic_Type=Simple
  Aggregation_Data_Set=None
Corner Coordinates:
Upper Left ( 0.0000000, 0.0000000)
Lower Left ( 0.000, 1800.000)
Upper Right ( 3600.000, 0.000)
Lower Right ( 3600.000, 1800.000)
Center ( 1800.000, 900.000)
Band 1 Block=64x64 Type=Int16, ColorInterp=Undefined
steph@bim:~/modis >

The problem seems to be with negative south and west values, as GRASS
needs south boundaries as 90S and not -90 (as my data have).

Is there a solution to use this kind of data in GRASS? Thanks for any
pointers on this topic.

Kind Regards

  Stephan Holl

--
Stephan Holl

Check headers for GnuPG Key!
http://www.gdf-hannover.de

2

On Sun, 14 Nov 2004, Stephan Holl wrote:

Dear Users,

I am trying to import some MODIS data into GRASS. In a first step I
managed to convert a subdataset of the MODIS data into img-format using
gdal_translate.

The dataset is a global dataset which should be in lat/long.

[...]

The problem seems to be with negative south and west values, as GRASS
needs south boundaries as 90S and not -90 (as my data have).

Is there a solution to use this kind of data in GRASS? Thanks for any
pointers on this topic.

What actually is the problem? Did something go wrong when you tried to import it into GRASS? Was this using r.in.gdal? And which GRASS version?

Paul

3

Hello Paul,

On Sun, 14 Nov 2004 11:55:34 +0000 (GMT) Paul Kelly
<paul-grass@stjohnspoint.co.uk> wrote:

On Sun, 14 Nov 2004, Stephan Holl wrote:

> Dear Users,
>
> I am trying to import some MODIS data into GRASS. In a first step I
> managed to convert a subdataset of the MODIS data into img-format
> using gdal_translate.
>
> The dataset is a global dataset which should be in lat/long.
>
[...]
>
>
> The problem seems to be with negative south and west values, as
> GRASS needs south boundaries as 90S and not -90 (as my data have).
>
> Is there a solution to use this kind of data in GRASS? Thanks for
> any pointers on this topic.

What actually is the problem? Did something go wrong when you tried to
import it into GRASS? Was this using r.in.gdal? And which GRASS
version?

yes, r.in.gdal did not work, because latitude for south was illigal.

Proceeding with import...
WARNING: G_set_window(): Illegal latitude for South

no map is imported.

thanks for any hints.

Stephan

--
Stephan Holl

Check headers for GnuPG Key!
http://www.gdf-hannover.de

4

When I was trying to expand my r.in.aster script to include MODIS data, I
got a message from Frank Warmerdam (GDAL head developer) that there was a
problem in georeferencing MODIS with GDAL. If I remember, I think it had to
do with the MODIS metadata that gdalwarp would use for georeferencing.

Michael Barton

On 11/14/04 11:31 AM, "Stephan Holl" <sholl@gmx.net> wrote:

Hello Paul,

On Sun, 14 Nov 2004 11:55:34 +0000 (GMT) Paul Kelly
<paul-grass@stjohnspoint.co.uk> wrote:

On Sun, 14 Nov 2004, Stephan Holl wrote:

Dear Users,

I am trying to import some MODIS data into GRASS. In a first step I
managed to convert a subdataset of the MODIS data into img-format
using gdal_translate.

The dataset is a global dataset which should be in lat/long.

[...]

The problem seems to be with negative south and west values, as
GRASS needs south boundaries as 90S and not -90 (as my data have).

Is there a solution to use this kind of data in GRASS? Thanks for
any pointers on this topic.

What actually is the problem? Did something go wrong when you tried to
import it into GRASS? Was this using r.in.gdal? And which GRASS
version?

yes, r.in.gdal did not work, because latitude for south was illigal.

Proceeding with import...
WARNING: G_set_window(): Illegal latitude for South

no map is imported.

thanks for any hints.

Stephan

____________________
C. Michael Barton, Professor of Anthropology
School of Human Evolution and Social Change
PO Box 872402
Arizona State University
Tempe, AZ 85287-2402
USA

Phone: 480-965-6262
Fax: 480-965-7671
www: <www.public.asu.edu/~cmbarton>

5

Hello Stefan

On Sun, 14 Nov 2004, Stephan Holl wrote:

What actually is the problem? Did something go wrong when you tried to
import it into GRASS? Was this using r.in.gdal? And which GRASS
version?

yes, r.in.gdal did not work, because latitude for south was illigal.

Proceeding with import...
WARNING: G_set_window(): Illegal latitude for South

no map is imported.

This doesn't really explain what looks like a bug, but I think it is a side-effect of GDAL trying to set the current region to match the imported data, which IMHO I don't think it should do. E.g. what if you were running an analysis in the background for which the region set was very important, and then imported a new file and r.in.gdal suddenly changed the current region?

I think it should only do this if the -e flag to extend the region is given, or perhaps there should be a new flag specifically to set the region to match the imported file.
If you comment out the lines in raster/r.in.gdal/main.c:
     if(G_set_window (&cellhd) < 0)
         exit(3);
and re-compile does it work?
Although seeing that doesn't fix the underlying bug, it might run into problems later I suppose.

Anybody else have an opinion on the region setting and extending behaviour of r.in.gdal and v.in.ogr?

Paul

6

Paul Kelly wrote:

>> What actually is the problem? Did something go wrong when you tried to
>> import it into GRASS? Was this using r.in.gdal? And which GRASS
>> version?
>
> yes, r.in.gdal did not work, because latitude for south was illigal.
>
> Proceeding with import...
> WARNING: G_set_window(): Illegal latitude for South
>
> no map is imported.

This doesn't really explain what looks like a bug, but I think it is a
side-effect of GDAL trying to set the current region to match the imported
data, which IMHO I don't think it should do. E.g. what if you were running
an analysis in the background for which the region set was very important,
and then imported a new file and r.in.gdal suddenly changed the current
region?

r.in.gdal doesn't modify the WIND file. If you use the -e flag, it
modifies the DEFAULT_WIND file in the PERMANENT mapset[1], otherwise
it only modifies the process' internal region.

[1] At least, it tries to; obviously, if you don't have write
permission of the PERMANENT mapset, this will fail.

I think it should only do this if the -e flag to extend the region is
given, or perhaps there should be a new flag specifically to set the
region to match the imported file.

I'm not sure that there would be much point, as you can just use
"g.region rast=..." if you want that behaviour.

If you comment out the lines in raster/r.in.gdal/main.c:
     if(G_set_window (&cellhd) < 0)
         exit(3);
and re-compile does it work?
Although seeing that doesn't fix the underlying bug, it might run into
problems later I suppose.

Anybody else have an opinion on the region setting and extending behaviour
of r.in.gdal and v.in.ogr?

If r.in.gdal is importing a georeferenced raster, it should attempt to
set the correct bounds and resolution (i.e. cellhd) for that raster,
rather than using (0,0,W,H). It shouldn't attempt to set the current
region or the default region; that should be left to g.region (which
should have an option to set the default region).

I don't think that this is an issue for v.in.ogr; AFAIK, vector maps
don't have a region (cellhd).

--
Glynn Clements <glynn@gclements.plus.com>

7

On Mon, 15 Nov 2004, Glynn Clements wrote:

Paul Kelly wrote:

This doesn't really explain what looks like a bug, but I think it is a
side-effect of GDAL trying to set the current region to match the imported
data, which IMHO I don't think it should do. E.g. what if you were running
an analysis in the background for which the region set was very important,
and then imported a new file and r.in.gdal suddenly changed the current
region?

r.in.gdal doesn't modify the WIND file. If you use the -e flag, it
modifies the DEFAULT_WIND file in the PERMANENT mapset[1], otherwise
it only modifies the process' internal region.

I don't disagree that this is the actual behaviour, but on reading though the source code it seems to me this wasn't the intended behaviour and it is trying to write out the bounds of the imported map as the current region (i.e WIND). But G_set_window() doesn't change anything because of some caching it does or something (I looked into it a long time ago and don't remember all the details now).

Unless I'm misunderstanding what G_set_window() is supposed to do but I think it's for setting the current region.

I think it should only do this if the -e flag to extend the region is
given, or perhaps there should be a new flag specifically to set the
region to match the imported file.

I'm not sure that there would be much point, as you can just use
"g.region rast=..." if you want that behaviour.

Yes I agree, except from what I've noticed it seems r.in.gdal was designed to actually do this. But maybe these lines should really be deleted then:

If you comment out the lines in raster/r.in.gdal/main.c:
     if(G_set_window (&cellhd) < 0)
         exit(3);
and re-compile does it work?
Although seeing that doesn't fix the underlying bug, it might run into
problems later I suppose.

Anybody else have an opinion on the region setting and extending behaviour
of r.in.gdal and v.in.ogr?

If r.in.gdal is importing a georeferenced raster, it should attempt to
set the correct bounds and resolution (i.e. cellhd) for that raster,
rather than using (0,0,W,H). It shouldn't attempt to set the current
region or the default region; that should be left to g.region (which
should have an option to set the default region).

I don't think that this is an issue for v.in.ogr; AFAIK, vector maps
don't have a region (cellhd).

If you are using v.in.ogr to create a new location based on the projection
of the vector map being imported, it has to create a DEFAULT_WIND, and it is
good to populate it with meaningful region extents, which can be got through
some OGR function (get extents or something like that).

Paul

8

Paul Kelly wrote:

>> This doesn't really explain what looks like a bug, but I think it is a
>> side-effect of GDAL trying to set the current region to match the imported
>> data, which IMHO I don't think it should do. E.g. what if you were running
>> an analysis in the background for which the region set was very important,
>> and then imported a new file and r.in.gdal suddenly changed the current
>> region?
>
> r.in.gdal doesn't modify the WIND file. If you use the -e flag, it
> modifies the DEFAULT_WIND file in the PERMANENT mapset[1], otherwise
> it only modifies the process' internal region.

I don't disagree that this is the actual behaviour, but on reading though
the source code it seems to me this wasn't the intended behaviour and it
is trying to write out the bounds of the imported map as the current
region (i.e WIND). But G_set_window() doesn't change anything because of
some caching it does or something (I looked into it a long time ago and
don't remember all the details now).

Unless I'm misunderstanding what G_set_window() is supposed to do but I
think it's for setting the current region.

G_set_window() sets the process' internal region (i.e. G__.window),
not the WIND file. It doesn't modify the WIND file, and isn't meant
to.

If you want to set the mapset's current region (i.e. the WIND file),
you need to use G_put_window(). If you want to modify some other file
(i.e. PERMANENT/DEFAULT_WIND or a cellhd file), you need to use
G__put_window().

Quite a few programs use G_set_window() to set a window other than
that defined in the WIND file.

E.g. programs which operate on the raw raster data without clipping,
padding or resampling (e.g. r.proj, r.bilinear) copy the raster's
cellhd to the current region:

  G_get_cellhd(mapname, mapset, &window);
  G_set_window(&window);

Similarly, import (r.in.*) programs set up a region with the same
number of rows and columns as the raster being imported. If
georeferencing information is available, they will set the boundaries
from that, otherwise they typically place the south-west corner at
(0,0) and the north-east corner at (W,H) (i.e. nsres == ewres == 1).

However, only those programs which wish to make persistent changes
(e.g. g.region, d.zoom) use G_put_window().

I'm fairly certain that r.in.gdal's behaviour (both with and without
the -e flag) is intentional.

>> I think it should only do this if the -e flag to extend the region is
>> given, or perhaps there should be a new flag specifically to set the
>> region to match the imported file.
>
> I'm not sure that there would be much point, as you can just use
> "g.region rast=..." if you want that behaviour.

Yes I agree, except from what I've noticed it seems r.in.gdal was designed
to actually do this. But maybe these lines should really be deleted then:

>> If you comment out the lines in raster/r.in.gdal/main.c:
>> if(G_set_window (&cellhd) < 0)
>> exit(3);

No. Any raster import program has to at least set the window such that
the rows/cols values match the actual number of rows and columns.
Either that, or it will have to manually resample the image.

G_put_raster_row() etc don't do any resampling. They expect the buffer
to contain <cols> values, and they expect to be called exactly <rows>
times.

Most raster programs take their region (and hence grid dimensions)
from the WIND file, but import programs take the grid dimensions from
the file being imported, and so have to set the region to match.

>> and re-compile does it work?
>> Although seeing that doesn't fix the underlying bug, it might run into
>> problems later I suppose.
>>
>> Anybody else have an opinion on the region setting and extending behaviour
>> of r.in.gdal and v.in.ogr?
>
> If r.in.gdal is importing a georeferenced raster, it should attempt to
> set the correct bounds and resolution (i.e. cellhd) for that raster,
> rather than using (0,0,W,H). It shouldn't attempt to set the current
> region or the default region; that should be left to g.region (which
> should have an option to set the default region).
>
> I don't think that this is an issue for v.in.ogr; AFAIK, vector maps
> don't have a region (cellhd).

If you are using v.in.ogr to create a new location based on the projection
of the vector map being imported, it has to create a DEFAULT_WIND, and it is
good to populate it with meaningful region extents, which can be got through
some OGR function (get extents or something like that).

Right; if you're creating a new location, it makes sense.

--
Glynn Clements <glynn@gclements.plus.com>

9

On Tue, 16 Nov 2004, Glynn Clements wrote:

Paul Kelly wrote:

[...]

I don't disagree that this is the actual behaviour, but on reading though
the source code it seems to me this wasn't the intended behaviour and it
is trying to write out the bounds of the imported map as the current
region (i.e WIND). But G_set_window() doesn't change anything because of
some caching it does or something (I looked into it a long time ago and
don't remember all the details now).

Unless I'm misunderstanding what G_set_window() is supposed to do but I
think it's for setting the current region.

G_set_window() sets the process' internal region (i.e. G__.window),
not the WIND file. It doesn't modify the WIND file, and isn't meant
to.

If you want to set the mapset's current region (i.e. the WIND file),
you need to use G_put_window(). If you want to modify some other file
(i.e. PERMANENT/DEFAULT_WIND or a cellhd file), you need to use
G__put_window().

Right; get you there.

[...]

I'm fairly certain that r.in.gdal's behaviour (both with and without
the -e flag) is intentional.

I wonder then why it doesn't just extend the current region, as (as I understand it) can be done reliably because the current user definitely has permission to change it, rather than extending the default.

If you are using v.in.ogr to create a new location based on the projection
of the vector map being imported, it has to create a DEFAULT_WIND, and it is
good to populate it with meaningful region extents, which can be got through
some OGR function (get extents or something like that).

Right; if you're creating a new location, it makes sense.

I don't think the -e flag is working in v.in.ogr. If I get a chance sometime I will try and update it to extend the current region, or perhaps it could be removed as regions aren't so relevant for vectors (should they be?)

Paul

10

Paul Kelly wrote:

[...]
> I'm fairly certain that r.in.gdal's behaviour (both with and without
> the -e flag) is intentional.

I wonder then why it doesn't just extend the current region, as (as I
understand it) can be done reliably because the current user definitely
has permission to change it, rather than extending the default.

Yeah, I thought that too.

Maybe it's because you can extend the current region yourself with
"g.region rast=...", but g.region doesn't allow you to change the
default region (which is a pretty major deficiency, IMHO).

--
Glynn Clements <glynn@gclements.plus.com>