Context Navigation

← Previous Ticket
Next Ticket →

#3193 closed enhancement (fixed)

r.out.gdal: add AUTHORITY node to the srs info of the exported raster

Reported by:	hellik	Owned by:	grass-dev@…
Priority:	major	Milestone:	7.6.0
Component:	Raster	Version:	svn-releasebranch76
Keywords:	r.external.out, r.out.gdal, v.external.out, v.out.ogr	Cc:
CPU:	All	Platform:	All

Description

test case in a location created by an EPSG code

g.proj -p                                                                       
-PROJ_INFO-------------------------------------------------
name       : MGI / Austria GK Central
datum      : hermannskogel
ellps      : bessel
proj       : tmerc
lat_0      : 0
lon_0      : 13.33333333333333
k          : 1
x_0        : 0
y_0        : -5000000
no_defs    : defined
towgs84    : 577.326,90.129,463.919,5.1366,1.4742,5.2970,2.4232
-PROJ_EPSG-------------------------------------------------
epsg       : 31255
-PROJ_UNITS------------------------------------------------
unit       : meter
units      : meters
meters     : 1

export a raster by r.out.gdal and do gdalinfo

gdalinfo raster_export_by_grass.tif
Driver: GTiff/GeoTIFF
Files: raster_export_by_grass.tif
Size is 7201, 7001
Coordinate System is:
PROJCS["MGI / Austria GK Central",
    GEOGCS["bessel",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel_1841",6377397.155,299.1528128],
            TOWGS84[577.326,90.129,463.919,5.1366,1.4742,5.297,2.4232]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",13.33333333333333],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",-5000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.3.svn with GDAL 2.1.1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -96005.000,  230005.000) ( 12d 3'57.90"E, 47d12' 8.96"N)
Lower Left  (  -96005.000,  159995.000) ( 12d 4'50.98"E, 46d34'22.11"N)
Upper Right (  -23995.000,  230005.000) ( 13d 0'59.64"E, 47d12'32.62"N)
Lower Right (  -23995.000,  159995.000) ( 13d 1'12.91"E, 46d34'45.25"N)
Center      (  -60000.000,  195000.000) ( 12d32'45.35"E, 46d53'30.77"N)

do gdal_translate and assign the srs

gdal_translate -a_srs EPSG:31255 raster_export_by_grass.tif geotif_by_gdaltranslate.tif

check by gdalinfo

gdalinfo geotif_by_gdaltranslate.tif
Driver: GTiff/GeoTIFF
Files: geotif_by_gdaltranslate.tif
       geotif_by_gdaltranslate.tif.aux.xml
Size is 7201, 7001
Coordinate System is:
PROJCS["MGI / Austria GK Central",
    GEOGCS["MGI",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232],
            AUTHORITY["EPSG","6312"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4312"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",13.33333333333333],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",-5000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","31255"]] <<==
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.3.svn with GDAL 2.1.1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -96005.000,  230005.000) ( 12d 3'57.90"E, 47d12' 8.96"N)
Lower Left  (  -96005.000,  159995.000) ( 12d 4'50.98"E, 46d34'22.11"N)
Upper Right (  -23995.000,  230005.000) ( 13d 0'59.64"E, 47d12'32.62"N)
Lower Right (  -23995.000,  159995.000) ( 13d 1'12.91"E, 46d34'45.25"N)
Center      (  -60000.000,  195000.000) ( 12d32'45.35"E, 46d53'30.77"N)

GDAL adds the AUTHORITY node of the corresponding EPSG code.

this may be help #3191 (and maybe also other GIS software reading GRASS exported data).

also mentioned in a GDAL ML thread:

From what I can see the driver doesn't handle adding a new entry in the 
spatial_ref_sys table if there's no AUTHORITY node, such as in your use case.
[...]
Make sure a EPSG code is attached to 
the input dataset, or assign it explictly

Change History (16)

comment:1 by sbl, 8 years ago

Keywords:	r.out.gdal r.external.out added

comment:2 by neteler, 8 years ago

Milestone:	7.2.0 → 7.2.1

Ticket retargeted after milestone closed

comment:3 by martinl, 8 years ago

Milestone:	7.2.1 → 7.2.2

comment:4 by martinl, 7 years ago

Milestone:	7.2.2 → 7.4.0

All enhancement tickets should be assigned to 7.4 milestone.

follow-up: 6 comment:5 by martinl, 7 years ago

This functionality was recently implemented in trunk, right? Reference is missing here.

in reply to: 5 ; follow-up: 9 comment:6 by mmetz, 7 years ago

Replying to martinl:

This functionality was recently implemented in trunk, right? Reference is missing here.

It was implemented by martinl with r71303

comment:7 by sbl, 7 years ago

Hmm... This is still an issue for me both in GRASS 7.4 (r72108) and GRASS 7.5 (r72157):
See also (sort of duplicate): #3048

g.proj -p
-PROJ_INFO-------------------------------------------------
name       : ETRS89 / UTM zone 33N
datum      : etrs89
ellps      : grs80
proj       : utm
zone       : 33
no_defs    : defined
-PROJ_EPSG-------------------------------------------------
epsg       : 25833
-PROJ_UNITS------------------------------------------------
unit       : meter
units      : meters
meters     : 1

r.external.out --verbose directory=/tmp extension=.tif format=GTiff options=COMPRESS=LZW
r.mapcalc expression="test=1"
gdalinfo /tmp/test.tif 

r.external.out -r
g.remove type=raster name=test -f
rm /tmp/test.tif
r.mapcalc expression="test=1"
r.out.gdal -c input=test output=/tmp/test.tif createopt="COMPRESS=LZW" --o
gdalinfo /tmp/test.tif

In both cases gdalinfo shows:

Coordinate System is:
PROJCS["UTM Zone 33, Northern Hemisphere",
    GEOGCS["grs80",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6258"]],
        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"]]]
Origin = (256495.000000000000000,6653695.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)

That means:

    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","25833"]]

is still missing at the bottom of the CRS information, compared to CRS info with:

gdal_edit.py -a_srs EPSG:25833 /tmp/test.tif
gdalinfo /tmp/test.tif

comment:8 by sbl, 7 years ago

This affects interoperability with e.g. GeoServer/GeoNode or QGIS (https://issues.qgis.org/issues/14856).

in reply to: 6 ; follow-up: 10 comment:9 by hellik, 7 years ago

Replying to mmetz:

Replying to martinl:

This functionality was recently implemented in trunk, right? Reference is missing here.

It was implemented by martinl with r71303

tested with the example above:

g.proj -p                                                                       
-PROJ_INFO-------------------------------------------------
name       : MGI / Austria GK Central
datum      : hermannskogel
ellps      : bessel
proj       : tmerc
lat_0      : 0
lon_0      : 13.33333333333333
k          : 1
x_0        : 0
y_0        : -5000000
no_defs    : defined
towgs84    : 577.326,90.129,463.919,5.1366,1.4742,5.2970,2.4232
-PROJ_EPSG-------------------------------------------------
epsg       : 31255
-PROJ_UNITS------------------------------------------------
unit       : meter
units      : meters
meters     : 1

r.out.gdal input=test1@test output=J:\test\test_write_auth.tif format=GTiff createopt=COMPRESS=LZW,PREDICTOR=2,TILED=YES

gdalinfo test_write_auth.tif
Driver: GTiff/GeoTIFF
Files: test_write_auth.tif
Size is 7201, 7001
Coordinate System is:
PROJCS["MGI / Austria GK Central",
    GEOGCS["bessel",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel_1841",6377397.155,299.1528128],
            TOWGS84[577.326,90.129,463.919,5.1366,1.4742,5.297,2.4232]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",13.33333333333333],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",-5000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.5.svn with GDAL 2.2.3
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -96005.000,  230005.000) ( 12d 3'57.90"E, 47d12' 8.96"N)
Lower Left  (  -96005.000,  159995.000) ( 12d 4'50.98"E, 46d34'22.11"N)
Upper Right (  -23995.000,  230005.000) ( 13d 0'59.64"E, 47d12'32.62"N)
Lower Right (  -23995.000,  159995.000) ( 13d 1'12.91"E, 46d34'45.25"N)
Center      (  -60000.000,  195000.000) ( 12d32'45.35"E, 46d53'30.77"N)

gdal_translate -a_srs EPSG:31255 test_write_auth.tif test_write_auth_assign_by_gdal.tif

Driver: GTiff/GeoTIFF
Files: test_write_auth_assign_by_gdal.tif
Size is 7201, 7001
Coordinate System is:
PROJCS["MGI / Austria GK Central",
    GEOGCS["MGI",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232],
            AUTHORITY["EPSG","6312"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4312"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",13.33333333333333],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",-5000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","31255"]] <==
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.5.svn with GDAL 2.2.3
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -96005.000,  230005.000) ( 12d 3'57.90"E, 47d12' 8.96"N)
Lower Left  (  -96005.000,  159995.000) ( 12d 4'50.98"E, 46d34'22.11"N)
Upper Right (  -23995.000,  230005.000) ( 13d 0'59.64"E, 47d12'32.62"N)
Lower Right (  -23995.000,  159995.000) ( 13d 1'12.91"E, 46d34'45.25"N)
Center      (  -60000.000,  195000.000) ( 12d32'45.35"E, 46d53'30.77"N)

System Info                                                                     
GRASS version: 7.5.svn                                                          
GRASS SVN revision: r72120                                                      
Build date: 2018-01-24                                                          
Build platform: i386-w64-mingw32                                                
GDAL: 2.2.3                                                                     
PROJ.4: 4.9.3                                                                   
GEOS: 3.5.0                                                                     
SQLite: 3.17.0                                                                  
Python: 2.7.4                                                                   
wxPython: 2.8.12.1                                                              
Platform: Windows-Vista-6.0.6002-SP2 (OSGeo4W)

comparing the gdalinfo of GRASS exported and GDAL assigned tif,

AUTHORITY["EPSG","31255"]

seems to be still missing

in reply to: 9 ; follow-ups: 12 14 comment:10 by mmetz, 7 years ago

Replying to hellik:

Replying to mmetz:

Replying to martinl:

This functionality was recently implemented in trunk, right? Reference is missing here.

It was implemented by martinl with r71303

tested with the example above: [...]
AUTHORITY["EPSG","31255"]
seems to be still missing

The functions were implemented in r71303. An EPSG code was only exported with v.out.ogr. As of trunk r72160, an existing EPSG code is exported with r.external.out, r.out.gdal, v.external.out, v.out.ogr.

comment:11 by neteler, 7 years ago

Milestone:	7.4.0 → 7.4.1

Ticket retargeted after milestone closed

in reply to: 10 ; follow-up: 13 comment:12 by hellik, 7 years ago

Replying to mmetz:

Replying to hellik:
Replying to mmetz:

Replying to martinl:

This functionality was recently implemented in trunk, right? Reference is missing here.

It was implemented by martinl with r71303

tested with the example above: [...]
AUTHORITY["EPSG","31255"]
seems to be still missing
The functions were implemented in r71303. An EPSG code was only exported with v.out.ogr. As of trunk r72160, an existing EPSG code is exported with r.external.out, r.out.gdal, v.external.out, v.out.ogr.

Backport candidat for 7.4.1?

in reply to: 12 comment:13 by mmetz, 7 years ago

Replying to hellik:

Replying to mmetz:
Replying to hellik:
Replying to mmetz:

Replying to martinl:

This functionality was recently implemented in trunk, right? Reference is missing here.

It was implemented by martinl with r71303

tested with the example above: [...]
AUTHORITY["EPSG","31255"]
seems to be still missing
The functions were implemented in r71303. An EPSG code was only exported with v.out.ogr. As of trunk r72160, an existing EPSG code is exported with r.external.out, r.out.gdal, v.external.out, v.out.ogr.
Backport candidat for 7.4.1?

It's a new feature, and EPSG codes should not be required in order to support custom coordinate reference systems.

in reply to: 10 comment:14 by hellik, 7 years ago

Replying to mmetz:

Replying to hellik:
Replying to mmetz:

Replying to martinl:

This functionality was recently implemented in trunk, right? Reference is missing here.

It was implemented by martinl with r71303

tested with the example above: [...]
AUTHORITY["EPSG","31255"]
seems to be still missing
The functions were implemented in r71303. An EPSG code was only exported with v.out.ogr. As of trunk r72160, an existing EPSG code is exported with r.external.out, r.out.gdal, v.external.out, v.out.ogr.

tested with

System Info                                                                     
GRASS version: 7.5.svn                                                          
GRASS SVN revision: r72194                                                      
Build date: 2018-01-30                                                          
Build platform: x86_64-w64-mingw32                                              
GDAL: 2.2.3                                                                     
PROJ.4: 4.9.3                                                                   
GEOS: 3.5.0                                                                     
SQLite: 3.17.0                                                                  
Python: 2.7.5                                                                   
wxPython: 2.8.12.1                                                              
Platform: Windows-8-6.2.9200 (OSGeo4W)

g.proj -p                                                                       
-PROJ_INFO-------------------------------------------------
name       : MGI / Austria GK Central
datum      : hermannskogel
ellps      : bessel
proj       : tmerc
lat_0      : 0
lon_0      : 13.33333333333333
k          : 1
x_0        : 0
y_0        : -5000000
no_defs    : defined
towgs84    : 577.326,90.129,463.919,5.1366,1.4742,5.2970,2.4232
-PROJ_EPSG-------------------------------------------------
epsg       : 31255
-PROJ_UNITS------------------------------------------------
unit       : meter
units      : meters
meters     : 1

gdalinfo test_raster_authority.tif
Driver: GTiff/GeoTIFF
Files: test_raster_authority.tif
Size is 170, 141
Coordinate System is:
PROJCS["MGI / Austria GK Central",
    GEOGCS["MGI",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[577.326,90.129,463.919,5.1366,1.4742,5.297,2.4232],
            AUTHORITY["EPSG","6312"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4312"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",13.33333333333333],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",-5000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","31255"]] <=
Origin = (-161500.000000000000000,273250.000000000000000)
Pixel Size = (250.000000000000000,-250.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.5.svn with GDAL 2.2.3
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -161500.000,  273250.000) ( 11d11'10.74"E, 47d34'42.28"N)
Lower Left  ( -161500.000,  238000.000) ( 11d11'56.99"E, 47d15'41.55"N)
Upper Right ( -119000.000,  273250.000) ( 11d45' 3.81"E, 47d35'15.35"N)
Lower Right ( -119000.000,  238000.000) ( 11d45'37.92"E, 47d16'14.25"N)
Center      ( -140250.000,  255625.000) ( 11d28'27.33"E, 47d25'29.61"N)

AUTHORITY["EPSG","XXXX"] is now exported in rasters too.

comment:15 by neteler, 7 years ago

Milestone:	7.4.1 → 7.4.2

comment:16 by neteler, 6 years ago

Keywords:	v.external.out v.out.ogr added
Milestone:	7.4.2 → 7.6.0
Resolution:	→ fixed
Status:	new → closed
Version:	svn-trunk → svn-releasebranch76

To be released in the upcoming 7.6.0; closing.

Note: See TracTickets for help on using tickets.

Download in other formats: