[GRASS-user] EPSG 3035

laurapoggio · March 25, 2009, 3:57pm

Dear all,
I am starting to use GRASS and I am trying to set up a location with the LAEA ETRS89 projection (EPSG 3035).
First I set up the location using the EPSG code and the I tried to import a geotiff with the same projection (created with gdal_translate). However I got the error “Projection of dataset does not appear to match current location”.
Then I set up the location using the georeferenced file itself, I tried to import the same file and I got the same error.

I can manage to import the file using the -o option and over-riding projection check, but I would like to understand and possible solve the problem.

I am using GRASS63 under linux.

Thank you very much in advance

Laura

marisn · March 25, 2009, 4:23pm

My guess - that GeoTIFF data are in EPSG:3035, still GeoTIFF's
coordinate reference system is wrong or missing - i.e. it has world
file and no proj file. GRASS will generate error in such case, still
if bounds are OK, importing with -o will do the right thing.

Laura, check Your GeoTIFF file with gdalinfo. If it says something like:
Coordinate System is `'
Origin = (722414.444181173457764,6363518.164592186920345)
then import into properly set location will work just fine.

Maris.

mlennert · March 25, 2009, 4:24pm

On 25/03/09 16:57, Laura Poggio wrote:

Dear all,
I am starting to use GRASS and I am trying to set up a location with the LAEA ETRS89 projection (EPSG 3035).
First I set up the location using the EPSG code and the I tried to import a geotiff with the same projection (created with gdal_translate). However I got the error "Projection of dataset does not appear to match current location".
Then I set up the location using the georeferenced file itself, I tried to import the same file and I got the same error.

I can manage to import the file using the -o option and over-riding projection check, but I would like to understand and possible solve the problem.

Your file might be georeferenced, but does not contain the information about the system it is georeferenced in. In that case, grass will always complain that it cannot guarantee a match.

If you are sure that it is referenced in 3035, then you have to use the -o option.

Moritz

laurapoggio · March 25, 2009, 5:05pm

The output of gdalinfo on the file is

Driver: GTiff/GeoTIFF
Files: data/JRC/map/DTM/ast1.tif
Size is 360, 364
Coordinate System is:
PROJCS[“ETRS89 / ETRS-LAEA”,
GEOGCS[“ETRS89”,
DATUM[“European_Terrestrial_Reference_System_1989”,
SPHEROID[“GRS 1980”,6378137,298.2572221010042,
AUTHORITY[“EPSG”,“7019”]],
AUTHORITY[“EPSG”,“6258”]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433],
AUTHORITY[“EPSG”,“4258”]],
UNIT[“metre”,1,
AUTHORITY[“EPSG”,“9001”]],
AUTHORITY[“EPSG”,“3035”]]
Origin = (4081404.118557849433273,2966229.602829793468118)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 4081404.119, 2966229.603)
Lower Left ( 4081404.119, 2955309.603)
Upper Right ( 4092204.119, 2966229.603)
Lower Right ( 4092204.119, 2955309.603)
Center ( 4086804.119, 2960769.603)
Band 1 Block=360x5 Type=Float32, ColorInterp=Gray
NoData Value=-32768
Metadata:
LAYER_TYPE=athematic

The parameters seems to be correct. I will then use the -o option

Thank you very much

Laura

2009/3/25 Moritz Lennert <mlennert@club.worldonline.be>

On 25/03/09 16:57, Laura Poggio wrote:

Dear all,
I am starting to use GRASS and I am trying to set up a location with the LAEA ETRS89 projection (EPSG 3035).
First I set up the location using the EPSG code and the I tried to import a geotiff with the same projection (created with gdal_translate). However I got the error “Projection of dataset does not appear to match current location”.
Then I set up the location using the georeferenced file itself, I tried to import the same file and I got the same error.

I can manage to import the file using the -o option and over-riding projection check, but I would like to understand and possible solve the problem.

Your file might be georeferenced, but does not contain the information about the system it is georeferenced in. In that case, grass will always complain that it cannot guarantee a match.

If you are sure that it is referenced in 3035, then you have to use the -o option.

Moritz

mlennert · March 25, 2009, 8:05pm

On 25/03/09 18:05, Laura Poggio wrote:

The output of gdalinfo on the file is

Driver: GTiff/GeoTIFF Files: data/JRC/map/DTM/ast1.tif Size is 360, 364 Coordinate System is: PROJCS["ETRS89 / ETRS-LAEA", GEOGCS["ETRS89", DATUM["European_Terrestrial_Reference_System_1989",
SPHEROID["GRS 1980",6378137,298.2572221010042,
AUTHORITY["EPSG","7019"]], AUTHORITY["EPSG","6258"]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","4258"]], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","3035"]]

Mmmh, so coordinate system info is included and GRASS should recognize it. Could you send the output of g.proj -w in the location you want to import the file into ?

Moritz

Paul_Kelly · March 25, 2009, 10:29pm

On Wed, 25 Mar 2009, Moritz Lennert wrote:

On 25/03/09 18:05, Laura Poggio wrote:

The output of gdalinfo on the file is

Driver: GTiff/GeoTIFF Files: data/JRC/map/DTM/ast1.tif Size is 360, 364 Coordinate System is: PROJCS["ETRS89 / ETRS-LAEA", GEOGCS["ETRS89", DATUM["European_Terrestrial_Reference_System_1989",
SPHEROID["GRS 1980",6378137,298.2572221010042,
AUTHORITY["EPSG","7019"]], AUTHORITY["EPSG","6258"]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","4258"]], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","3035"]]

Mmmh, so coordinate system info is included and GRASS should recognize it.

The WKT is invalid though as it says it's projected, but doesn't include either the name of the projection nor any parameters for it.

Paul

neteler · March 26, 2009, 6:54am

On Wed, Mar 25, 2009 at 11:29 PM, Paul Kelly
<paul-grass@stjohnspoint.co.uk> wrote:

On Wed, 25 Mar 2009, Moritz Lennert wrote:

On 25/03/09 18:05, Laura Poggio wrote:

The output of gdalinfo on the file is

Driver: GTiff/GeoTIFF Files:
data/JRC/map/DTM/ast1.tif Size is 360, 364 Coordinate
System is: PROJCS["ETRS89 / ETRS-LAEA",
GEOGCS["ETRS89",
DATUM["European_Terrestrial_Reference_System_1989",
SPHEROID["GRS 1980",6378137,298.2572221010042,
AUTHORITY["EPSG","7019"]], AUTHORITY["EPSG","6258"]],
PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4258"]], UNIT["metre",1,
AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","3035"]]

Mmmh, so coordinate system info is included and GRASS should recognize it.

The WKT is invalid though as it says it's projected, but doesn't include
either the name of the projection nor any parameters for it.

But if I used OGR's tool, I can generate it, too:

PYTHONPATH=/usr/local/lib64/python2.5/site-packages/osgeo epsg_tr.py 3035
PROJCS["ETRS89 / ETRS-LAEA",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",52],
    PARAMETER["longitude_of_center",10],
    PARAMETER["false_easting",4321000],
    PARAMETER["false_northing",3210000],
    AUTHORITY["EPSG","3035"],
    AXIS["X",EAST],
    AXIS["Y",NORTH]]

grep 'ETRS89 / ETRS-LAEA' lib/proj/pcs.csv
3035,"ETRS89 / ETRS-LAEA",9001,4258,19986,9820,1,0,8802,10,9102,8801,52,9102,8807,3210000,9001,8806,4321000,9001,

?
Markus

laurapoggio · March 26, 2009, 7:39am

The output of the g.proj -w command is:

WARNING: Unable to initialise PROJ.4 with the following parameter list:
+a=6378137 +rf=298.257222101 +no_defs +towgs84=0.000,0.000,0.000
WARNING: The error message: projection not named
WARNING: Can’t parse GRASS PROJ_INFO file
WARNING: g.proj: Unable to convert to WKT

Then I think there is some problem

Thank you very much

Laura

2009/3/25 Moritz Lennert <mlennert@club.worldonline.be>

On 25/03/09 18:05, Laura Poggio wrote:

The output of gdalinfo on the file is

Driver: GTiff/GeoTIFF Files: data/JRC/map/DTM/ast1.tif Size is 360, 364 Coordinate System is: PROJCS[“ETRS89 / ETRS-LAEA”, GEOGCS[“ETRS89”, DATUM[“European_Terrestrial_Reference_System_1989”,
SPHEROID[“GRS 1980”,6378137,298.2572221010042,
AUTHORITY[“EPSG”,“7019”]], AUTHORITY[“EPSG”,“6258”]], PRIMEM[“Greenwich”,0], UNIT[“degree”,0.0174532925199433], AUTHORITY[“EPSG”,“4258”]], UNIT[“metre”,1, AUTHORITY[“EPSG”,“9001”]], AUTHORITY[“EPSG”,“3035”]]

Mmmh, so coordinate system info is included and GRASS should recognize it. Could you send the output of g.proj -w in the location you want to import the file into ?

Moritz

mlennert · March 26, 2009, 8:32am

On 26/03/09 08:39, Laura Poggio wrote:

The output of the g.proj -w command is:

WARNING: Unable to initialise PROJ.4 with the following parameter list:
+a=6378137 +rf=298.257222101 +no_defs +towgs84=0.000,0.000,0.000
WARNING: The error message: projection not named
WARNING: Can't parse GRASS PROJ_INFO file
WARNING: g.proj: Unable to convert to WKT

Then I think there is some problem

What about g.proj -p ?

Moritz

mlennert · March 26, 2009, 9:24am

On 26/03/09 10:03, Laura Poggio wrote:

The output of g.proj - p

-PROJ_INFO-------------------------------------------------
datum : etrs89
ellps : grs80
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1

That looks better

No, you should have something like this:

GRASS 7.0.svn (test3035):~ > g.proj -p
-PROJ_INFO-------------------------------------------------
name : Lambert Azimuthal Equal Area
proj : laea
datum : etrs89
ellps : grs80
lat_0 : 52
lon_0 : 10
x_0 : 4321000
y_0 : 3210000
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1
GRASS 7.0.svn (test3035):~ > g.proj -w
PROJCS["Lambert Azimuthal Equal Area",
     GEOGCS["grs80",
         DATUM["European_Terrestrial_Reference_System_1989",
SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]],
     PROJECTION["Lambert_Azimuthal_Equal_Area"],
     PARAMETER["latitude_of_center",52],
     PARAMETER["longitude_of_center",10],
     PARAMETER["false_easting",4321000],
     PARAMETER["false_northing",3210000],
     UNIT["metre",1]]

How did you create this location ?

Moritz

mlennert · March 26, 2009, 9:41am

[please keep the discussion on the mailing list]

On 26/03/09 10:34, Laura Poggio wrote:

I used the "create from georeferenced file" and selected the file I posted the gdalinfo in this mail.

Ok, as Paul pointed out the projection information in that file is insufficient, so it is logical that GRASS cannot create the correct location, but it should at least say so during the process and probably refuse to create the location.

I have the same output also from the location created using EPSG code and the value 3035.

This is strange, you should get below results. Which versions of GRASS, proj and gdal are you using ?

Moritz

Thank you

Laura

2009/3/26 Moritz Lennert <mlennert@club.worldonline.be <mailto:mlennert@club.worldonline.be>>

    On 26/03/09 10:03, Laura Poggio wrote:

        The output of g.proj - p

        -PROJ_INFO-------------------------------------------------
        datum : etrs89
        ellps : grs80
        -PROJ_UNITS------------------------------------------------
        unit : metre
        units : metres
        meters : 1

        That looks better

    No, you should have something like this:

    GRASS 7.0.svn (test3035):~ > g.proj -p
    -PROJ_INFO-------------------------------------------------
    name : Lambert Azimuthal Equal Area
    proj : laea
    datum : etrs89
    ellps : grs80
    lat_0 : 52
    lon_0 : 10
    x_0 : 4321000
    y_0 : 3210000
    no_defs : defined

    -PROJ_UNITS------------------------------------------------
    unit : metre
    units : metres
    meters : 1
    GRASS 7.0.svn (test3035):~ > g.proj -w
    PROJCS["Lambert Azimuthal Equal Area",
       GEOGCS["grs80",

           DATUM["European_Terrestrial_Reference_System_1989",

    SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
           PRIMEM["Greenwich",0],
           UNIT["degree",0.0174532925199433]],

       PROJECTION["Lambert_Azimuthal_Equal_Area"],
       PARAMETER["latitude_of_center",52],
       PARAMETER["longitude_of_center",10],
       PARAMETER["false_easting",4321000],
       PARAMETER["false_northing",3210000],
       UNIT["metre",1]]

    How did you create this location ?

    Moritz

laurapoggio · March 26, 2009, 10:04am

Sorry I hit reply instead of reply to all.

The versions of the software are:
GRASS 6.3.0
GDAL 1.5.3
PROJ 4.5.0
These are the latest available on the repositories of fedora9. But I will have a further look and try to update at least GDAL.

I deleted the existing locations and did everything from the beginning in order to be able to copy the exact messages.

Creation of location with EPSG code:
Datum not recognized by GRASS and no parameters found.

Creation of location with georeferenced file (the one I posted details):
no error/warning messages

Thank you

laura

2009/3/26 Moritz Lennert <mlennert@club.worldonline.be>

[please keep the discussion on the mailing list]

On 26/03/09 10:34, Laura Poggio wrote:

I used the “create from georeferenced file” and selected the file I posted the gdalinfo in this mail.

Ok, as Paul pointed out the projection information in that file is insufficient, so it is logical that GRASS cannot create the correct location, but it should at least say so during the process and probably refuse to create the location.

I have the same output also from the location created using EPSG code and the value 3035.

This is strange, you should get below results. Which versions of GRASS, proj and gdal are you using ?

Moritz

Thank you

Laura

2009/3/26 Moritz Lennert <mlennert@club.worldonline.be mailto:[mlennert@club.worldonline.be](mailto:mlennert@club.worldonline.be)>

On 26/03/09 10:03, Laura Poggio wrote:

The output of g.proj - p

-PROJ_INFO-------------------------------------------------
datum : etrs89
ellps : grs80
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1

That looks better

No, you should have something like this:

GRASS 7.0.svn (test3035):~ > g.proj -p
-PROJ_INFO-------------------------------------------------
name : Lambert Azimuthal Equal Area
proj : laea
datum : etrs89
ellps : grs80
lat_0 : 52
lon_0 : 10
x_0 : 4321000
y_0 : 3210000
no_defs : defined

-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1
GRASS 7.0.svn (test3035):~ > g.proj -w
PROJCS[“Lambert Azimuthal Equal Area”,
GEOGCS[“grs80”,

DATUM[“European_Terrestrial_Reference_System_1989”,

SPHEROID[“Geodetic_Reference_System_1980”,6378137,298.257222101]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433]],

PROJECTION[“Lambert_Azimuthal_Equal_Area”],
PARAMETER[“latitude_of_center”,52],
PARAMETER[“longitude_of_center”,10],
PARAMETER[“false_easting”,4321000],
PARAMETER[“false_northing”,3210000],
UNIT[“metre”,1]]

How did you create this location ?

Moritz

mlennert · March 27, 2009, 9:52am

On 26/03/09 11:04, Laura Poggio wrote:

Sorry I hit reply instead of reply to all.

The versions of the software are:
GRASS 6.3.0

Ok, trying with 6.3 and proj 4.6, I get:

GRASS 6.3.0 (newLocation):~ > g.proj -c location=eee epsg=3035
WARNING: Datum <unknown> not recognised by GRASS and no parameters found.
Location eee created!

And g.proj -p in that new location gives:

-PROJ_INFO-------------------------------------------------
name : Lambert Azimuthal Equal Area
proj : laea
ellps : wgs84
lat_0 : 52
lon_0 : 10
x_0 : 4321000
y_0 : 3210000
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : Meter
units : Meters
meters : 1

However, trying the same in GRASS 6.5svn and same proj version, I get no warning and g.proj -p in the resulting location gives:

g.proj -p
-PROJ_INFO-------------------------------------------------
name : Lambert Azimuthal Equal Area
proj : laea
datum : etrs89
ellps : grs80
lat_0 : 52
lon_0 : 10
x_0 : 4321000
y_0 : 3210000
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1

So, don't know exactly where the problem lies but it seems to be fixed in later versions (can't try 6.4rc right now). Cc'ing Paul in case he knows.

Moritz

laurapoggio · March 27, 2009, 12:05pm

I will then try to upgrade the version of GRASS, but in the meantime just a question.

All the files I will work on have the same projection as the location, at least according to gdal and using EPSG codes.
If I create the location with the problems discussed until now and then import the files with -o option, would this be source of troubles for further elaborations such as watershed analysis?

Thank you very much

Laura

2009/3/27 Moritz Lennert <mlennert@club.worldonline.be>

On 26/03/09 11:04, Laura Poggio wrote:

Sorry I hit reply instead of reply to all.

The versions of the software are:
GRASS 6.3.0

Ok, trying with 6.3 and proj 4.6, I get:

GRASS 6.3.0 (newLocation):~ > g.proj -c location=eee epsg=3035
WARNING: Datum not recognised by GRASS and no parameters found.
Location eee created!

And g.proj -p in that new location gives:

-PROJ_INFO-------------------------------------------------
name : Lambert Azimuthal Equal Area
proj : laea

ellps : wgs84
lat_0 : 52
lon_0 : 10
x_0 : 4321000
y_0 : 3210000
no_defs : defined
-PROJ_UNITS------------------------------------------------

unit : Meter
units : Meters
meters : 1

However, trying the same in GRASS 6.5svn and same proj version, I get no warning and g.proj -p in the resulting location gives:

g.proj -p
-PROJ_INFO-------------------------------------------------
name : Lambert Azimuthal Equal Area
proj : laea
datum : etrs89
ellps : grs80
lat_0 : 52
lon_0 : 10
x_0 : 4321000
y_0 : 3210000
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1

So, don’t know exactly where the problem lies but it seems to be fixed in later versions (can’t try 6.4rc right now). Cc’ing Paul in case he knows.

Moritz

mlennert · March 27, 2009, 12:29pm

On 27/03/09 13:05, Laura Poggio wrote:

I will then try to upgrade the version of GRASS, but in the meantime just a question.

All the files I will work on have the same projection as the location, at least according to gdal and using EPSG codes.
If I create the location with the problems discussed until now and then import the files with -o option, would this be source of troubles for further elaborations such as watershed analysis?

I don't think so. But you might want to try to create a "correct" location with all the parameters, including datum, etc. You might have to do this by hand though, i.e. using the option to create the location with projection parameters, instead of EPSG code or existing file.

Moritz

laurapoggio · March 27, 2009, 12:56pm

Ok thanks. Maybe by hand is really the easiest solution.

Thank you

Laura

2009/3/27 Moritz Lennert <mlennert@club.worldonline.be>

On 27/03/09 13:05, Laura Poggio wrote:

I will then try to upgrade the version of GRASS, but in the meantime just a question.

All the files I will work on have the same projection as the location, at least according to gdal and using EPSG codes.
If I create the location with the problems discussed until now and then import the files with -o option, would this be source of troubles for further elaborations such as watershed analysis?

I don’t think so. But you might want to try to create a “correct” location with all the parameters, including datum, etc. You might have to do this by hand though, i.e. using the option to create the location with projection parameters, instead of EPSG code or existing file.

Moritz