r.out.gdal does not write EPSG code of the projection

[sbl]

GeoTiffs produced with GRASS GIS 6.4, 7.0.4 and 7.1.svn do not contain the AUTHORITY parameter for the projection, which is thus not properly recognized by other tools (e.g. QGIS, but esp. GeoServer).

Related discussion on the GDAL mailinglist can be found here: ​http://thread.gmane.org/gmane.comp.gis.gdal.devel/43019

g.proj -w gave the projection information as in the GeoTIFF, though without any AUTHORITY parameters (as it is supposed to according to the manual). The AUTHORITY tags seem to be added by r.out.gdal, except for the authority parameter for the entire projection. Even if g.proj -g returns epsg=25832.

Compared to GDALs GeoTIFFs also AXIS parameters are missing (see below).

GRASS 7.0.4svn (ETRS_32N):~ >g.proj -g
name=ETRS89 / UTM zone 32N
datum=etrs89
ellps=grs80
proj=utm
zone=32
towgs84=0,0,0,0,0,0,0
no_defs=defined
epsg=25832
unit=meter
units=meters
meters=1
GRASS 7.0.4svn (ETRS_32N):~ >
GRASS 7.0.4svn (ETRS_32N):~ >r.out.gdal input=test output=$HOME/test.tif
GRASS 7.0.4svn (ETRS_32N):~ >gdalinfo $HOME/test.tif
(...)
Coordinate System is:
PROJCS["UTM Zone 32, Northern Hemisphere",
    GEOGCS["grs80",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (0.000000000000000,1.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
(...)
GRASS 7.0.4svn (ETRS_32N):~ >
GRASS 7.0.4svn (ETRS_32N):~ >gdal_translate -a_srs EPSG:25832 $HOME/test.tif $HOME/test2.tif
GRASS 7.0.4svn (ETRS_32N):~ >gdalinfo $HOME/test2.tif
(...)
Coordinate System is:
PROJCS["ETRS89 / UTM zone 32N",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","25832"]]
Origin = (0.000000000000000,1.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
(...)

[neteler]

I just tried with the NC dataset on my system, no such issue:

GRASS 7.0.5svn (nc_spm_08_grass7):~ > r.out.gdal elevation output=elevation.tif
Checking GDAL data type and nodata value...
 100%
Using GDAL data type <Float32>
Exporting raster data to GTiff format...
ERROR 6: SetColorTable() only supported for Byte or UInt16 bands in TIFF format.
 100%
r.out.gdal complete. File <elevation.tif> created.
GRASS 7.0.5svn (nc_spm_08_grass7):~ > gdalinfo elevation.tif 
Driver: GTiff/GeoTIFF
Files: elevation.tif
Size is 1500, 1350
Coordinate System is:
PROJCS["Lambert Conformal Conic",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",36.16666666666666],
    PARAMETER["standard_parallel_2",34.33333333333334],
    PARAMETER["latitude_of_origin",33.75],
    PARAMETER["central_meridian",-79],
    PARAMETER["false_easting",609601.22],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (630000.000000000000000,228500.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  630000.000,  228500.000) ( 78d46'27.41"W, 35d48'34.59"N)
Lower Left  (  630000.000,  215000.000) ( 78d46'28.63"W, 35d41'16.54"N)
Upper Right (  645000.000,  228500.000) ( 78d36'29.89"W, 35d48'33.08"N)
Lower Right (  645000.000,  215000.000) ( 78d36'32.01"W, 35d41'15.03"N)
Center      (  637500.000,  221750.000) ( 78d41'29.49"W, 35d44'54.91"N)
Band 1 Block=1500x1 Type=Float32, ColorInterp=Gray
  Metadata:
    COLOR_TABLE_RULES_COUNT=5
    COLOR_TABLE_RULE_RGB_0=5.557880e+01 7.572900e+01 0 191 191 0 255 0
    COLOR_TABLE_RULE_RGB_1=7.572900e+01 9.587920e+01 0 255 0 255 255 0
    COLOR_TABLE_RULE_RGB_2=9.587920e+01 1.160290e+02 255 255 0 255 127 0
    COLOR_TABLE_RULE_RGB_3=1.160290e+02 1.361800e+02 255 127 0 191 127 63
    COLOR_TABLE_RULE_RGB_4=1.361800e+02 1.563300e+02 191 127 63 200 200 200
    Generated_with=GRASS GIS 7.0.5svn

GRASS 7.0.5svn (nc_spm_08_grass7):~ > gdal-config --version
2.0.2

But doing the same in a UTM32N location, the AUTHORITY is missing:

GRASS 7.0.5svn (utm32n):~ > g.proj -w
PROJCS["UTM Zone 32, Northern Hemisphere",
    GEOGCS["wgs84",
        DATUM["WGS_1984",
            SPHEROID["WGS_1984",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

# pure GDAL:
epsg_tr.py 32632
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32632"]]

I see differences in the PROJCS, perhaps the text "warping" partially fails?

[neteler]

Thinking more about it, I believe that it is not possible to generate the AUTHORITY here because GRASS internally uses the PROJ.4 representation.

Options:

• implement lookup support in GDAL, see ​http://article.gmane.org/gmane.comp.gis.gdal.devel/43031
• Change libproj to first read PROJ_EPSG (if present), then merge the datum information into that which is maintained in GRASS only, send all to GDAL and get WKT back (something like this).

Note that EPSG codes often don't include the DATUM string since several datums may exist (that's why GRASS comes with a sophisticated lookup mechanism which isn't present in GDAL for example).

Back to the original report:

which is thus not properly recognized by other tools (e.g. QGIS, but esp. GeoServer).

Can you elaborate?

[sbl]

To me it seems the issue is also present in your nc_spm_08_grass7 location...

Replying to neteler:

Thinking more about it, I believe that it is not possible to generate the AUTHORITY here because GRASS internally uses the PROJ.4 representation.

In my location the AUTHORITY parameter is present but it is not used in r.out.gdal never the less, so it would not have to be generated.

Note that EPSG codes often don't include the DATUM string since several datums may exist (that's why GRASS comes with a sophisticated lookup mechanism which isn't present in GDAL for example).

Back to the original report:

which is thus not properly recognized by other tools (e.g. QGIS, but esp. GeoServer).

Can you elaborate?

Sure, in QGIS my raster in EPSG:25832 is interpreted as EPSG:3044 (both have equal proj-strings). However, CRS do not match and also in PostGIS SRID would not match, so they have to defined manually for analysis with other data one has in the originally used CRS 25832. In GeoNode / GeoServer the GeoTiffs exported from GRASS are shown as pink boxes as long as one does not define the projection in GeoServer manually. See: ​http://docs.geonode.org/en/master/tutorials/advanced/geonode_production/adv_gsconfig/crs_handling.html

[sbl]

See also possibly related tickets for QGIS processing: ​https://hub.qgis.org/issues/11884 and ​https://hub.qgis.org/issues/14499

[sbl]

This is also true for r.external.out, see: #3059

[mmetz]

Replying to sbl:

This is also true for r.external.out, see: #3059

Granted that the file PERMANENT/PROJ_EPSG exists, this works for me with relbr74 and trunk, no action required.

Explanation:

The EPSG code of the projection is only written out if it is known. Use g.proj -p or g.proj -g to test if the EPSG code is known.

Regarding GDAL/OGR input: if the EPSG code of the projection of the input is known and if a new location is created from that input, this EPSG code is now preserved in PERMANENT/PROJ_EPSG as of trunk r72114.

[sbl]

Closing as duplicate. See: #3193