[GRASS-dev] GDAL write support

[Glynn Clements]

Markus Neteler 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...

Done in r41521.

> 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


If I use:

	r.external.out
	r.mapcalc 'foo = float(elevation.dem)'

The resulting map works within GRASS:

	Driver: GTiff/GeoTIFF
	Files: /opt/grass-data/spearfish57/glynn/gdal/foo
	Size is 1899, 1398
	Coordinate System is `'
	Origin = (590010.000000000000000,4928000.000000000000000)
	Pixel Size = (10.000000000000000,-10.000000000000000)
	Image Structure Metadata:
	  INTERLEAVE=BAND
	Corner Coordinates:
	Upper Left  (  590010.000, 4928000.000) 
	Lower Left  (  590010.000, 4914020.000) 
	Upper Right (  609000.000, 4928000.000) 
	Lower Right (  609000.000, 4914020.000) 
	Center      (  599505.000, 4921010.000) 
	Band 1 Block=1899x1 Type=Float32, ColorInterp=Gray
	  NoData Value=nan

> > 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?

It's decided at open time.

For integer maps, it's GDT_Byte, GDT_UInt16 or GDT_Int32, as set by
G_set_cell_format(), defaulting to GDT_Int32.

For FP maps, it's GDT_Float32 for FCELL and GDT_Float64 for DCELL (the
choice is determined by Rast_set_fp_type(), defaulting to DCELL if
$GRASS_FP_DOUBLE exists and FCELL otherwise).


> > 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"));

The code already prints a message if it uses the memory driver;
lib/raster/gdal.c, lines 448-455:

    /* If not - create MEM driver for intermediate dataset. 
     * Check if raster can be created at all (with GDALCreateCopy) */
    else if ((*pGDALGetMetadataItem) (driver, GDAL_DCAP_CREATECOPY, NULL)) {
	GDALDriverH mem_driver;

	G_message(_("Driver <%s> does not support direct writing. "
		    "Using MEM driver for intermediate dataset."),
		  st->opts.format);

> > 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?

It applies to anything; I don't know how r.out.gdal itself handles it.

> > 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.

I'm not sure that you understood what I wrote. The code to convert
GRASS' projection information to WKT cannot be used by lib/raster,
only by modules, so it cannot be used at the point that a raster map
is written via GDAL.

A WKT version could be created by r.external.out then simply copied
into the output at creation time. But if you change the projection,
the WKT data will then be out of sync. If you want to be able to write
maps with projection information via GDAL, g.proj really needs to be
changed to maintain a WKT version of the projection information as
well as the GRASS version.

-- 
Glynn Clements <glynn at gclements.plus.com>