[GRASS-user] Loosing projection info?

[Hamish]

iomeneandrei wrote:
> I have a location with these default region settings:
> 
> projection: 1 (UTM)
> zone:       33
> datum:      eur50
> ellipsoid:  international
....
> When I export this raster with r.out.gdal input=dem at debug format=GTiff
> 'output=C:/tif.tif', I have a geotiff file that
> seems to have lost some projection related metainfo.
> 
> If I apply gdalinfo to this exported tif I have:
> [...]
> Driver: GTiff/GeoTIFF
> Size is 765, 552
> Coordinate System is:
> PROJCS["UTM Zone 33, Northern Hemisphere",
>     GEOGCS["international",
>         DATUM["unknown",
>            
> SPHEROID["unnamed",6378388,297.0000000000014]],
>         PRIMEM["Greenwich",0],
>         UNIT["degree",0.0174532925199433]],
>     PROJECTION["Transverse_Mercator"],
>     PARAMETER["latitude_of_origin",0],
>     PARAMETER["central_meridian",15],
>     PARAMETER["scale_factor",0.9996],
>     PARAMETER["false_easting",500000],
>     PARAMETER["false_northing",0],
>     UNIT["metre",1,
>         AUTHORITY["EPSG","9001"]]]
> [...]
> 
> In Grass I have a defined datum, but it's undefined in the output
> geotiff.

what does "g.proj -p", "g.proj -w", and "g.proj -j" say for you?


> For the same reason probably I have a wrong epsg code (the
> right one should be 23033).

epsg 9001 has to do with the definition of the meter, not your map.
(follow the [,[,]] pairs)

 
> Is it normal, or is there something wrong in my procedure?

the problem is not a mistake by you AFAICT--
 (but g.proj output please, and what version of grass and gdal?)

for OGR there is a known problem with datum grid files, perhaps the
issue exists with GDAL rasters as well?   see this thread from last week:
  http://thread.gmane.org/gmane.comp.gis.grass.devel/30047

but AFAIK that was just for +nadgrid files, other 3,7 term +towgs84
params should make it into the GeoTIFF.


hmmm, I just created a new location based on EPSG 23033 using the GRASS
startup GUI, selected the generic ("all of eu50 zone", second on list)
datum option, and exported a dummy raster with r.out.gdal.

gdalinfo shows:

Coordinate System is:
PROJCS["UTM Zone 33, Northern Hemisphere",
    GEOGCS["international",
        DATUM["European_Datum_1950",
            SPHEROID["International 1924",6378388,297.000000000005,
                AUTHORITY["EPSG","7022"]],
            AUTHORITY["EPSG","6230"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]


ie datum and ellipsoid are now named, but +towgs84 params are missing.


GRASS64> g.proj -p
-PROJ_INFO-------------------------------------------------
name       : Universe Transverse Mercator
proj       : utm
datum      : eur50
ellps      : international
zone       : 33
no_defs    : defined
towgs84    : -87.000,-98.000,-121.000
-PROJ_UNITS------------------------------------------------
unit       : metre
units      : metres
meters     : 1

GRASS64> g.proj -j
+proj=utm
+no_defs
+zone=33
+a=6378388
+rf=297
+towgs84=-87.000,-98.000,-121.000
+to_meter=1

GRASS64> g.proj -w
PROJCS["UTM Zone 33, Northern Hemisphere",
    GEOGCS["international",
        DATUM["European_Datum_1950",
            SPHEROID["International_1924",6378388,297],
            TOWGS84[-87.000,-98.000,-121.000]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1]]

note that here TOWGS84 *is* listed in the WKT.



Hamish