[GRASS5] GRASS 5.3/5.7, GDAL and southern hemisphere

GRASS GRASS developer
,
1

Today we imported a LANDSAT scene into a southern
hemisphere location (Argentina), UTM. There is a slight
confusion with 'south: defined' and negative zone numbers:

gdalinfo tells us:

GRASS 5.7.0-cvs:~ > gdalinfo tuc6.tif
Driver: GTiff/GeoTIFF
Size is 1024, 692
Coordinate System is:
PROJCS["WGS 84 / UTM zone 4S",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235629972,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-159],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32704"]]
Origin = (-65.356600,9999973.275500)
Pixel Size = (30.000000,-30.000000)
Corner Coordinates:
Upper Left ( -65.357, 9999973.275)
Lower Left ( -65.357, 9979213.275)
Upper Right ( 30654.643, 9999973.275)
Lower Right ( 30654.643, 9979213.275)
Center ( 15294.643, 9989593.275)
Band 1 Block=1024x1 Type=Byte, ColorInterp=Gray

Importing it:

GRASS 5.7.0-cvs:~ > r.in.gdal in=tuc6.tif out=tuc6
100%
CREATING SUPPORT FILES FOR tuc6
SETTING GREY COLOR TABLE FOR tuc6 (8bit, full range)

Error on display:

GRASS 5.7.0-cvs:~ > d.rast tuc6
WARNING: [tuc6] in mapset [tucuman] - in different zone [-4] than current
         region [4]
WARNING: unable to open raster map [tuc6 in tucuman]
ERROR: Not able to open cellfile for [tuc6]

How the location was defined:

GRASS 5.7.0-cvs:~ > g.projinfo -p

-----------------------------------------------------------
PROJ_INFO file:
name: UTM
datum: wgs84
towgs84: 0.000,0.000,0.000
proj: utm
ellps: wgs84
a: 6378137.0000000000
es: 0.0066943800
f: 298.2572235630
zone: 4
south: defined
-----------------------------------------------------------
PROJ_UNITS file:
unit: meter
units: meters
meters: 1.0
-----------------------------------------------------------

I am not sure where it should be fixed.
Thanks for hints,

Markus

2

Hello Markus

On Thu, 25 Sep 2003, Markus Neteler wrote:

Today we imported a LANDSAT scene into a southern
hemisphere location (Argentina), UTM. There is a slight
confusion with 'south: defined' and negative zone numbers:

gdalinfo tells us:

GRASS 5.7.0-cvs:~ > gdalinfo tuc6.tif
Driver: GTiff/GeoTIFF
Size is 1024, 692
Coordinate System is:
PROJCS["WGS 84 / UTM zone 4S",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235629972,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-159],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32704"]]
Origin = (-65.356600,9999973.275500)
Pixel Size = (30.000000,-30.000000)
Corner Coordinates:
Upper Left ( -65.357, 9999973.275)
Lower Left ( -65.357, 9979213.275)
Upper Right ( 30654.643, 9999973.275)
Lower Right ( 30654.643, 9979213.275)
Center ( 15294.643, 9989593.275)
Band 1 Block=1024x1 Type=Byte, ColorInterp=Gray

Importing it:

GRASS 5.7.0-cvs:~ > r.in.gdal in=tuc6.tif out=tuc6
100%
CREATING SUPPORT FILES FOR tuc6
SETTING GREY COLOR TABLE FOR tuc6 (8bit, full range)

Error on display:

GRASS 5.7.0-cvs:~ > d.rast tuc6
WARNING: [tuc6] in mapset [tucuman] - in different zone [-4] than current
         region [4]
WARNING: unable to open raster map [tuc6 in tucuman]
ERROR: Not able to open cellfile for [tuc6]

Seems like maybe the location is not set up properly? Did you do it with
g.setproj or r.in.gdal. Could you try editing the DEFAULT_WIND and WIND
files manually to change the zone to -4? Does that help?

3

On Thu, 25 Sep 2003, Markus Neteler wrote:

WARNING: [tuc6] in mapset [tucuman] - in different zone [-4] than current
         region [4]

Looking into it further, it appears r.in.gdal will set the zone in the
raster header to a negative value if the region is in the southern
hemisphere. But g.setproj will never create a negative zone in the WIND
file (which has the same format as a raster/CELL header), and indicates
that it is in the southern hemisphere by putting south: defined in the
PROJ_INFO file.

pj_get_kv() will accept either south: defined or a negative zone number in
the PROJ_INFO. I suspect if you created the location using r.in.gdal it
would put a negative zone number in the WIND, and d.rast would be
happy.

To summarise: this inconsistency could be fixed either in r.in.gdal or
g.setproj and I'm not sure which would be more appropriate as I don't
understand the significancy of being in the southern hemisphere for the
UTM projection. Maybe somebody could explain this.

Temporary workaround: make sure the zone in the cellhd and WIND are the
same---I'm not sure if they're actually used for anything these days
anyway.

Paul

4

Hello Paul,

thanks for your answers.

On Thu, Sep 25, 2003 at 03:26:53PM +0100, Paul Kelly wrote:

On Thu, 25 Sep 2003, Markus Neteler wrote:
> WARNING: [tuc6] in mapset [tucuman] - in different zone [-4] than current
> region [4]

So - it appears that the colleague has TM and not UTM (so partially
solved). But:

Looking into it further, it appears r.in.gdal will set the zone in the
raster header to a negative value if the region is in the southern
hemisphere. But g.setproj will never create a negative zone in the WIND
file (which has the same format as a raster/CELL header), and indicates
that it is in the southern hemisphere by putting south: defined in the
PROJ_INFO file.

pj_get_kv() will accept either south: defined or a negative zone number in
the PROJ_INFO. I suspect if you created the location using r.in.gdal it
would put a negative zone number in the WIND, and d.rast would be
happy.

I assume that it is always happy when the map zone number matches
the overall defined zone number.

To summarise: this inconsistency could be fixed either in r.in.gdal or
g.setproj and I'm not sure which would be more appropriate as I don't
understand the significancy of being in the southern hemisphere for the
UTM projection. Maybe somebody could explain this.

Maybe Hamish can help?

Temporary workaround: make sure the zone in the cellhd and WIND are the
same---I'm not sure if they're actually used for anything these days
anyway.

Will do, thanks.

Markus

5

> To summarise: this inconsistency could be fixed either in r.in.gdal
> or g.setproj and I'm not sure which would be more appropriate as I
> don't understand the significancy of being in the southern
> hemisphere for the UTM projection. Maybe somebody could explain
> this.

Maybe Hamish can help?

The difference is that UTM quadrants in the Southern Hemisphere have a
false northing of 10,000,000. This is to prevent the possibility of
negative map coordinates.

Hamish

p.s.-
If it is any use, the following Matlab (or gpl'd Octave) script
will tell you what UTM zone you are in for a given [lat lon].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% utm_zone.m Calc UTM zone from Lat/Lon position.
%% Hamish Bowman, mid 2002 sometime
%
% UTM is based on a Transverse Mercator (conformal, cylindrical) projection.
% 60 zones. Width 6deg lon, numbered 1 to 60, starting at 180deg lon (West).
% 8deg lat strips, C to X northwards, omitting I and O, beginning at 80deg
% South.
%
% Panama: Clark1866, NAD27
% Zone:

% West is negative, south is negative
if (exist('LatLon') ~= 1)
        LatLon = [7.871289533 -81.6107378 ]
end

Strips = [ 'C' 'D' 'E' 'F' 'G' 'H' 'J' 'K' 'L' 'M' ...
                   'N' 'P' 'Q' 'R' 'S' 'T' 'U' 'V' 'W' 'X'];

Zone = ceil((LatLon(2)+180)/6)
Strip = Strips(ceil((LatLon(1)+80)/8))

6

On Thu, 25 Sep 2003, Paul Kelly wrote:

On Thu, 25 Sep 2003, Markus Neteler wrote:

> WARNING: [tuc6] in mapset [tucuman] - in different zone [-4] than current
> region [4]

Looking into it further, it appears r.in.gdal will set the zone in the
raster header to a negative value if the region is in the southern
hemisphere. But g.setproj will never create a negative zone in the WIND
file (which has the same format as a raster/CELL header), and indicates
that it is in the southern hemisphere by putting south: defined in the
PROJ_INFO file.

And looking even further, it appears that some old coorcnv library
functions add a large false northing if the zone in the WIND file is
negative, so it appears this was the original intention of the GRASS
authors and the authors of g.setproj (which was not developed at CERL if I
remember correctly) just got this bit wrong.

So I have changed g.setproj to put a negative zone number in the
DEFAULT_WIND if the answer is yes to the southern hemisphere question. But
probably in the future this issue is bound to cause a few more problems to
people because of the inconsistencies we have observed between different
GRASS modules.

Paul

7

On Tue, Oct 07, 2003 at 11:48:28AM +0100, Paul Kelly wrote:

On Thu, 25 Sep 2003, Paul Kelly wrote:

>
>
> On Thu, 25 Sep 2003, Markus Neteler wrote:
>
> > WARNING: [tuc6] in mapset [tucuman] - in different zone [-4] than current
> > region [4]
>
> Looking into it further, it appears r.in.gdal will set the zone in the
> raster header to a negative value if the region is in the southern
> hemisphere. But g.setproj will never create a negative zone in the WIND
> file (which has the same format as a raster/CELL header), and indicates
> that it is in the southern hemisphere by putting south: defined in the
> PROJ_INFO file.

And looking even further, it appears that some old coorcnv library
functions add a large false northing if the zone in the WIND file is
negative, so it appears this was the original intention of the GRASS
authors and the authors of g.setproj (which was not developed at CERL if I
remember correctly) just got this bit wrong.

So I have changed g.setproj to put a negative zone number in the
DEFAULT_WIND if the answer is yes to the southern hemisphere question. But
probably in the future this issue is bound to cause a few more problems to
people because of the inconsistencies we have observed between different
GRASS modules.

Paul,

thanks for the update. It seems to work OK now (Andrey Kinselev was
so kind to add Transv. Mercator support to LANDSAT/FAST driver in
GDAL which is needed for Argentina data).

So it seems to be complete now.

Markus