r.external.out does not write CRS information

[sbl]

At least for GeoTiffs, which is the format I used, files created with r.external.out do not contain any CRS information.

The CRS is a central piece of information for a spatial data set and should not get lost during export.

In order to reproduce try:

g.proj -g
r.external.out directory=tmp extension=.tif format=GTiff options="COMPRESS=LZW,PREDICTOR=2"
r.mapcalc expression="test=1"
gdalinfo $GISDBASE/$LOCATION/$MAPSET/tmp/test.tif | grep -A 5 "Coordinate System"

[sbl]

Some additional observations: A colleague and me are facing significant "interoperability problems" with GeoTiffs produced with r.external.out, which (after FTP transfer) could not be opened in ENVI.

ENVI complains about broken corner coordinates and a corrupted LZW-table.

The files are solar radiation grids (all integer) and are saved as Int32 by r.external.out. They were results from a tiled processing which was merged using r.patch. When I export via r.out.gdal I get a warning message that Metadata exceeds 32000 bytes and cannot be written to GeoTiff, and therefore will be transferred to PAM. So this is possibly related to #3056

[neteler]

I can reproduce it with the NC location:

# NC Location
g.region raster=elevation -p

g.proj -w
PROJCS["Lambert Conformal Conic",
    GEOGCS["grs80",
        DATUM["North_American_Datum_1983",
            SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    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["Meter",1]]

mkdir -p $HOME/tmp
r.external.out directory=$HOME/tmp extension=.tif format=GTiff options="COMPRESS=LZW,PREDICTOR=2"

r.mapcalc expression="test = 1"

gdalinfo $HOME/tmp/test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 1500, 1350
Coordinate System is `'
Origin = (630000.000000000000000,228500.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
...

Apparently some GDAL part is not initialized completely. I see in

​https://trac.osgeo.org/grass/browser/grass/trunk/lib/raster/gdal.c#L431

that the use of GPJ_grass_to_wkt() is commented. Maybe that part needs to be updated/completed?

[glynn]

Replying to neteler:

Apparently some GDAL part is not initialized completely. I see in

​https://trac.osgeo.org/grass/browser/grass/trunk/lib/raster/gdal.c#L431

that the use of GPJ_grass_to_wkt() is commented. Maybe that part needs to be updated/completed?

The main problem is that gproj depends upon libgis, which means that libgis cannot depend upon gproj. However, it may be possible to load gproj dynamically, in the same manner as GDAL. The only other option is to merge gproj into libgis, making GDAL (along with its large and open-ended set of dependencies) a hard dependency of libgis.

[neteler]

Replying to glynn:

Replying to neteler:

Apparently some GDAL part is not initialized completely. I see in

​https://trac.osgeo.org/grass/browser/grass/trunk/lib/raster/gdal.c#L431

that the use of GPJ_grass_to_wkt() is commented. Maybe that part needs to be updated/completed?

The main problem is that gproj depends upon libgis, which means that libgis cannot depend upon gproj. However, it may be possible to load gproj dynamically, in the same manner as GDAL.

Rather than merging libproj into libgis, your suggestion of dynamical loading of libproj sounds more viable to avoid a dependency "hell".

[sbl]

See also: #3191 and #3193

[neteler]

Since we run again into this problem, is there a way to "attach" the CRS?

[mmetz]

Replying to neteler:

Since we run again into this problem, is there a way to "attach" the CRS?

Yes, there is no dependency hell any more because library raster fns are now (as of GRASS 7) in a separate library which can nowbe dependent on libproj.

Please try trunk r72098 and relbr74 r72099.

[sbl]

Thanks Markus! I can confirm that projection information now is written to GeoTiffs produced with r.external.out in GRASS 7.4. However, like in r.out.gdal EPSG codes of the projections are not written here either (see #3048). Even for Locations purely created from EPSG code. Would be really helpful if EPSG codes somehow could be added where appropriate (and recorded in the projection information in GRASS).

[mmetz]

adjust milestone