[GRASS-dev] GDAL write support

[Markus Neteler]

On Sat, Dec 6, 2008 at 10:27 AM, Glynn Clements
<glynn at gclements.plus.com> wrote:
>
> I have committed preliminary, minimal, experimental, barely-tested
> support for creating new maps via GDAL (i.e. like r.external, but in
> the opposite direction).
>
> The feature is enabled with r.external.out, which allows you to
> optionally set the format (default GeoTIFF), extension (default none;
> use e.g. extension=.tif to set it), directory where output files are
> stored (default: <mapset>/gdal/) and any creation options (passed
> directly to the driver).
>
> Thereafter, any new raster maps will be generated via GDAL. The maps
> should automatically be available to GRASS, as if they had been linked
> with r.external.

Wish: could a flag the added to print the current state? Say, I am
activating it in a mapset, then stop working and reuse the mapset
after some month, I hardly remember how things were...

I have tested the module today:

 r.external terra_lst1km20030314.LST_Day_1km.tif out=modis_celsius
 # define output directory for GRASS calculation results:
 r.external.out $HOME/gisoutput/

 # do something (here: extract pixels > 20°C), write output directly as GeoTIFF:
 r.mapcalc "warm.tif = if(modis_celsius > 20, modis_celsius, null() )"

 # use the result elsewhere
 qgis $HOME/gisoutput/warm.tif

The resulting map, unfortunately is all NULL:

gdalinfo -mm $HOME/gisoutput/warm.tif
...
Band 1 Block=3264x1 Type=Int32, ColorInterp=Gray
    Computed Min/Max=-2147483648.000,-2147483648.000
  NoData Value=2147483648


It works ok with integer:

 r.mapcalc "myint.tif = col()"
gdalinfo -mm myint.tif
...
 Band 1 Block=34x27 Type=Int32, ColorInterp=Gray
    Computed Min/Max=1.000,34.000
  NoData Value=2147483648


Another FLOAT test:

r.mapcalc "myfloat.tif = 0.5 / col() "
qgis myfloat.tif
-> partially broken as nodata strips appear


> Running "r.external.out -r" will revert to creating maps natively.
> Maps previously created via GDAL should still work.
>
> Limitations include (but probably aren't limited to):
>
> 1. The data type for integer maps is determined by
> G_set_cell_format(), meaning that it will normally be Int32. Unlike
> r.out.gdal, we don't know the range of the data when the output file
> is opened.

As I understand it, we cannot know the map type a priori. Is it changed
on the fly in case of float maps?

> 2. If the driver for the format doesn't support direct writing, the
> data will intially be written via the memory driver, then copied to
> the file upon close. I'm guessing that this requires the map to fit
> into memory. r.external.out should print a warning if you choose such
> a driver.

Is this test reasonable?

    if( GDALGetDriverByName( "MEM" ) != NULL )
        G_warning(_("MEMORY driver is used which can lead to swapping"));


> 3. If you write data type which the format can't handle (e.g. writing
> Float64 to PNG), you lose.

... this should also apply to r.out.gdal, right?

> 4. Colour tables aren't included in the file. GRASS treats colour
> tables as distinct entities; the code which creates the raster data
> doesn't know about the colour table, or if one will even be created.
> Colour tables are normally only written after the map is closed, at
> which point it may be too late to add it to the image file.
>
> 5. Projection information isn't written either. The code which would
> need to do this is in lib/gis, which can't use e.g. GPJ_grass_to_wkt()
> from lib/proj because of circular dependency issues. I suppose that I
> could make r.external.out create a WKT version of the projection data,
> but then you run into problems if the projection is subsequently
> changed.

IMHO a WKT file would be better than losing the information completely.

> 6. You can't store multiple maps in a single multi-band file.

Not a severe problem since GDAL tools are available for that conversion.

The module will be very attractive to those who want to use GRASS as
GIS backbone while keeping the data outside of GRASS in common
formats.

Markus