r.external - was: Re: [GRASS-dev] some questions about future development

GRASS GRASS developer
,
1

Glynn,

this is excellent! Thanks so much. I have tested it but no success yet:

On Wed, Aug 20, 2008 at 1:46 PM, Glynn Clements
<glynn@gclements.plus.com> wrote:

Glynn Clements wrote:
> > > BTW, a r.external would be a great addition.

It's not that much work.

...

Tomorrow, I'll look into creating an r.external module, adding null
support, etc.

Okay, done (well, enough to try it out without having to manually hack
the files; null support is absent).

[for my convenience I have copied it over to GRASS 6.4.svn locally]

GRASS 6.4.svn (pat):~ > gdalinfo
/geodata2_originals_raid5/pat_OFD_RGB_IT2006_GB/060100.tif
Driver: GTiff/GeoTIFF
Files: /geodata2_originals_raid5/pat_OFD_RGB_IT2006_GB/060100.tif
Size is 13440, 11760
Coordinate System is:
PROJCS["Monte Mario / Italy zone 1",
    GEOGCS["Monte Mario",
        DATUM["Monte_Mario",
            SPHEROID["International 1924",6378388,297.0000000000014,
                AUTHORITY["EPSG","7022"]],
            AUTHORITY["EPSG","6265"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4265"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","3003"]]
Origin = (1660920.000000000000000,5107440.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=Adobe Photoshop CS2 Macintosh
  TIFFTAG_DATETIME=2007:03:01 06:12:15
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 1660920.000, 5107440.000) ( 11d 4'54.85"E, 46d 6'2.53"N)
Lower Left ( 1660920.000, 5101560.000) ( 11d 4'47.69"E, 46d 2'52.14"N)
Upper Right ( 1667640.000, 5107440.000) ( 11d10'7.61"E, 46d 5'56.71"N)
Lower Right ( 1667640.000, 5101560.000) ( 11d10'0.15"E, 46d 2'46.33"N)
Center ( 1664280.000, 5104500.000) ( 11d 7'27.57"E, 46d 4'24.46"N)
Band 1 Block=13440x1 Type=Byte, ColorInterp=Red
  Overviews: 6720x5880, 3360x2940, 1680x1470, 840x735, 420x368, 210x184
Band 2 Block=13440x1 Type=Byte, ColorInterp=Green
  Overviews: 6720x5880, 3360x2940, 1680x1470, 840x735, 420x368, 210x184
Band 3 Block=13440x1 Type=Byte, ColorInterp=Blue
  Overviews: 6720x5880, 3360x2940, 1680x1470, 840x735, 420x368, 210x184

GRASS 6.4.svn (pat):~ > r.external
/geodata2_originals_raid5/pat_OFD_RGB_IT2006_GB/060100.tif
out=ortho_trento
Projection of input dataset and current location appear to match
<ortho_trento> created
r.external complete.

GRASS 6.4.svn (pat):~ > g.region rast=ortho_trento -p
projection: 99 (Transverse Mercator)
zone: 0
datum: rome40
ellipsoid: international
north: 5107440
south: 5101560
west: 1660920
east: 1667640
nsres: 0.5
ewres: 0.5
rows: 11760
cols: 13440
cells: 158054400

GRASS 6.4.svn (pat):~ > d.rast ortho_trento

WARNING: Error reading map <ortho_trento@PERMANENT>, row 0
100%

r.info -r ortho_trento
WARNING: category support for [ortho_trento] in mapset [PERMANENT] missing
min=NULL
max=NULL

GRASS 6.4.svn (pat):~ > r.support ortho_trento
Edit header for [ortho_trento]? (y/n) [n]
Update the statistics (histogram, range) for [ortho_trento]? (y/n) [n] y

Updating statistics for [ortho_trento]

WARNING: Error reading map <ortho_trento@PERMANENT>, row 0
WARNING: Histogram for [ortho_trento in PERMANENT] missing (run r.support)
Edit the category file for [ortho_trento]? (y/n) [n]
Create/Update the color table for [ortho_trento]? (y/n) [n]
Edit the history file for [ortho_trento]? (y/n) [n]

The null file for [ortho_trento] may indicate that some cells contain
no data. If the null file for [ortho_trento] doesn't exist, zero cells in
it are treated by GRASS application programs as no data.

Do you want to create/reset the null file for [ortho_trento] so that
null cell values are considered valid data? (y/n) [n]

Do you want to delete the null file for [ortho_trento]
(all zero cells will then be considered no data)? (y/n) [n]

GRASS 6.4.svn (pat):~ > ls -la grassdata/pat/PERMANENT/cell_misc/ortho_trento/
total 86
drwxrwxr-x 2 neteler neteler 1024 2008-08-20 13:49 .
drwx---r-x 8 neteler neteler 84992 2008-08-20 13:49 ..
-rw-rw-r-- 1 neteler neteler 92 2008-08-20 13:49 gdal
-rw-rw-r-- 1 neteler neteler 0 2008-08-20 13:49 range

GRASS 6.4.svn (pat):~ > cat grassdata/pat/PERMANENT/cell_misc/ortho_trento/gdal
file: /geodata2_originals_raid5/pat_OFD_RGB_IT2006_GB/060100.tif
band: 1
null: none
type: 1

GRASS 6.4.svn (pat):~ > cat grassdata/pat/PERMANENT/cellhd/ortho_trento
proj: 99
zone: 0
north: 5107440
south: 5101560
east: 1667640
west: 1660920
cols: 13440
rows: 11760
e-w resol: 0.5
n-s resol: 0.5
format: 0
compressed: 0

The same happens with a GeoTIFF without overviews..
And the same happens in GRASS 7.

Markus

2

Markus Neteler wrote:

this is excellent! Thanks so much. I have tested it but no success yet:

>> Tomorrow, I'll look into creating an r.external module, adding null
>> support, etc.
>
> Okay, done (well, enough to try it out without having to manually hack
> the files; null support is absent).

[for my convenience I have copied it over to GRASS 6.4.svn locally]

Copied what over? r.external won't work unless the corresponding
lib/gis changes have also been incorporated, and lib/gis was built
with GDAL_LINK=1:

At present, r.external is built automatically, but the code to read
linked maps is conditionalised upon GDAL_LINK, so you need to build
with "make GDAL_LINK=1". If you forget that part, r.external will
work, but attempting to read the map will fail (the cell/fcell files
are zero-length).

Your errors are consistent with lib/gis not having the necessary
support. Does:

  nm $GISBASE/lib/libgrass_gis.so | fgrep gdal

produce any output? There should be at least:

  T G_close_gdal_link
  T G_get_gdal_link
  t gdal_values_double
  t gdal_values_float
  t gdal_values_int
  t read_data_gdal

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

3

On Wed, Aug 20, 2008 at 3:46 PM, Glynn Clements
<glynn@gclements.plus.com> wrote:

Markus Neteler wrote:

this is excellent! Thanks so much. I have tested it but no success yet:

>> Tomorrow, I'll look into creating an r.external module, adding null
>> support, etc.
>
> Okay, done (well, enough to try it out without having to manually hack
> the files; null support is absent).

[for my convenience I have copied it over to GRASS 6.4.svn locally]

Copied what over? r.external won't work unless the corresponding
lib/gis changes have also been incorporated, and lib/gis was built
with GDAL_LINK=1:

At present, r.external is built automatically, but the code to read
linked maps is conditionalised upon GDAL_LINK, so you need to build
with "make GDAL_LINK=1". If you forget that part, r.external will
work, but attempting to read the map will fail (the cell/fcell files
are zero-length).

Bingo. I overlooked the GDAL_LINK thing.

Your errors are consistent with lib/gis not having the necessary
support.

Right.

Now r.stats works! Will try to (locally) backport to 6.4 now (for
monitor tests).

Thanks again,
Markus

PS: Small issue:

gcc -I/home/neteler/grass70/dist.x86_64-unknown-linux-gnu/include -g
-Wall -fno-common -mtune=nocona -m64 -minline-all-stringops -fPIC
-D_FILE_OFFSET_BITS=64 -DGDAL_LINK -DPACKAGE=\""grasslibs"\"
-I/usr/local/include -I/usr/local/include
-I/home/neteler/grass70/dist.x86_64-unknown-linux-gnu/include -o
OBJ.x86_64-unknown-linux-gnu/quant.o -c quant.c
quant.c:293:1: warning: "MIN" redefined
In file included from /usr/local/include/gdal.h:40,
                 from
/home/neteler/grass70/dist.x86_64-unknown-linux-gnu/include/grass/gis.h:32,
                 from quant.c:286:
/usr/local/include/cpl_port.h:257:1: warning: this is the location of
the previous definition
quant.c:294:1: warning: "MAX" redefined
/usr/local/include/cpl_port.h:258:1: warning: this is the location of
the previous definition

4

Markus Neteler wrote:

>> this is excellent! Thanks so much. I have tested it but no success yet:
>
>> >> Tomorrow, I'll look into creating an r.external module, adding null
>> >> support, etc.
>> >
>> > Okay, done (well, enough to try it out without having to manually hack
>> > the files; null support is absent).
>>
>> [for my convenience I have copied it over to GRASS 6.4.svn locally]
>
> Copied what over? r.external won't work unless the corresponding
> lib/gis changes have also been incorporated, and lib/gis was built
> with GDAL_LINK=1:
>
>> At present, r.external is built automatically, but the code to read
>> linked maps is conditionalised upon GDAL_LINK, so you need to build
>> with "make GDAL_LINK=1". If you forget that part, r.external will
>> work, but attempting to read the map will fail (the cell/fcell files
>> are zero-length).

Bingo. I overlooked the GDAL_LINK thing.

> Your errors are consistent with lib/gis not having the necessary
> support.

Right.

Now r.stats works! Will try to (locally) backport to 6.4 now (for
monitor tests).

Thanks again,
Markus

PS: Small issue:
quant.c:293:1: warning: "MIN" redefined
In file included from /usr/local/include/gdal.h:40,
                 from
/home/neteler/grass70/dist.x86_64-unknown-linux-gnu/include/grass/gis.h:32,
                 from quant.c:286:
/usr/local/include/cpl_port.h:257:1: warning: this is the location of
the previous definition
quant.c:294:1: warning: "MAX" redefined
/usr/local/include/cpl_port.h:258:1: warning: this is the location of
the previous definition

Fixed, also the one in fpreclass.c.

I #undef'd any existing definition, as the GDAL version is broken, as
were the versions in quant.c, fpreclass.c, r.external, and r.in.gdal.

Now that I mention it, also see broken versions of MAX() or max() in:

  lib/external/shapelib/shpopen.c
  lib/g3d/tilewrite.c
  raster/r.carve/enforce_ds.c
  raster/wildfire/r.ros/spot_dist.c (*)
  raster3d/base/r3.mask.main.c
  raster3d/base/r3.null.main.c
  raster3d/r3.in.ascii/main.c
  raster3d/r3.in.v5d/main.c
  raster3d/r3.out.ascii/main.c
  raster3d/r3.out.v5d/main.c
  vector/v.in.ogr/main.c

(*) r.ros is broken in the opposite way to the rest.

When defining macros, *always* put parentheses around the expression
and all of the the arguments. IOW:

  #define max(a,b) ((a) < (b) ? (a) : (b))

not (like in r.ros):

  #define max(a,b) (a) < (b) ? (a) : (b)

or (like all of the others):

  #define max(a,b) (a < b ? a : b)

Otherwise, it will do the wrong thing if the arguments are complex
expressions, or the macro is used within a complex expression. E.g.:

  #define max(a,b) (a < b ? a : b)

  max(x & 0xFF, y & 0xFF)
== x & 0xFF < y & 0xFF ? x & 0xFF : y & 0xFF
== x & (0xFF < y) & 0xFF ? x & 0xFF : y & 0xFF

[as < has higher precedence than &]

Or:

  #define max(a,b) (a) < (b) ? (a) : (b)

  max(x, y) + 1
== (x) < (y) ? (x) : (y) + 1
== (x) < (y) ? (x) : ((y) + 1)

[as + has higher precedence than ?:]

Even if it doesn't actually matter when you write the macro (e.g.
because the arguments are variables and the macro call is the entire
RHS of an assignment), you never know when someone will add a "+ 1" or
whatever and spend hours figuring out what's wrong. Macros may look
like functions, but they don't behave anything like them.

Also, h2g2 reference:

  #include <stdio.h>
  #define SIX 1 + 5
  #define NINE 8 + 1
  int main(void) {
      printf("SIX * NINE = %d\n", SIX * NINE);
      return 0;
  }
=> 42

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

5

On Wed, Aug 20, 2008 at 5:34 PM, Markus Neteler <neteler@osgeo.org> wrote:

On Wed, Aug 20, 2008 at 3:46 PM, Glynn Clements
<glynn@gclements.plus.com> wrote:

Markus Neteler wrote:

...

Copied what over? r.external won't work unless the corresponding
lib/gis changes have also been incorporated, and lib/gis was built
with GDAL_LINK=1:

I have now backported all to GRASS 6.4.svn (locally) including the fixes
of the last 2 hours (MAX etc fixes).

...

Now r.stats works! Will try to (locally) backport to 6.4 now (for
monitor tests).

Done.
COOL STUFF. I linked a GeoTIFF RGB orthophoto (600MB)
and d.rgb works sufficiently fast.

Any objections that I submit the backport?

I assume that later the GDAL_LINK condition will become default?

Markus

6

Markus Neteler wrote:

>> Copied what over? r.external won't work unless the corresponding
>> lib/gis changes have also been incorporated, and lib/gis was built
>> with GDAL_LINK=1:

I have now backported all to GRASS 6.4.svn (locally) including the fixes
of the last 2 hours (MAX etc fixes).

...
> Now r.stats works! Will try to (locally) backport to 6.4 now (for
> monitor tests).

Done.
COOL STUFF. I linked a GeoTIFF RGB orthophoto (600MB)
and d.rgb works sufficiently fast.

Any objections that I submit the backport?

I assume that later the GDAL_LINK condition will become default?

It's reasonable enough to make it the default, although I really would
like to support a --without-gdal option (specifically, the ability to
only use GDAL for r.{in,out}.gdal and v.{in,out}.ogr, so that if GDAL
doesn't work, you only lose those modules, not the whole of GRASS).

AFAICT, the only "core" GDAL dependencies are the OSR stuff in
lib/proj and g.proj, and the vector library's support for v.external.
In which case, it really ought to be possible to get a usable version
of GRASS with no GDAL dependency in any of the libraries or core
modules (e.g. g.proj, g.region).

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

7

On Wed, 20 Aug 2008, Glynn Clements wrote:

AFAICT, the only "core" GDAL dependencies are the OSR stuff in
lib/proj and g.proj, and the vector library's support for v.external.
In which case, it really ought to be possible to get a usable version
of GRASS with no GDAL dependency in any of the libraries or core
modules (e.g. g.proj, g.region).

It used to be possible to compile g.proj without GDAL support, and use it only for reporting on the projection information and verifying datum information, but I removed all the conditional stuff for that when GDAL became a non-optional dependency. Could be put back in I suppose.

Paul

8

Paul Kelly wrote:

> AFAICT, the only "core" GDAL dependencies are the OSR stuff in
> lib/proj and g.proj, and the vector library's support for v.external.
> In which case, it really ought to be possible to get a usable version
> of GRASS with no GDAL dependency in any of the libraries or core
> modules (e.g. g.proj, g.region).

It used to be possible to compile g.proj without GDAL support, and use it
only for reporting on the projection information and verifying datum
information, but I removed all the conditional stuff for that when GDAL
became a non-optional dependency. Could be put back in I suppose.

That would be useful.

Also, is there any fundamental obstacle to being able to set the
projection information using PROJ parameters without going via OSR?
Given that GRASS uses PROJ for the actual projection, it shouldn't
actually need anything other than PROJ parameters. g.setproj doesn't
use any lib/proj functions except GPJ_ask_datum_params() (which no
longer exists).

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

9

On Thu, 21 Aug 2008, Glynn Clements wrote:

Paul Kelly wrote:

AFAICT, the only "core" GDAL dependencies are the OSR stuff in
lib/proj and g.proj, and the vector library's support for v.external.
In which case, it really ought to be possible to get a usable version
of GRASS with no GDAL dependency in any of the libraries or core
modules (e.g. g.proj, g.region).

It used to be possible to compile g.proj without GDAL support, and use it
only for reporting on the projection information and verifying datum
information, but I removed all the conditional stuff for that when GDAL
became a non-optional dependency. Could be put back in I suppose.

That would be useful.

I'll get on to it but it won't be before the end of the weekend I'd say.

Also, is there any fundamental obstacle to being able to set the
projection information using PROJ parameters without going via OSR?
Given that GRASS uses PROJ for the actual projection, it shouldn't
actually need anything other than PROJ parameters. g.setproj doesn't
use any lib/proj functions except GPJ_ask_datum_params() (which no
longer exists).

The main issue is the short ellipsoid names, which can be quite different. I believe the GRASS ellipse.table pre-dates the integration of PROJ into GRASS in the early 1990s. OSR understands the PROJ ellipsoid names, so at present it is used to get the underlying ellipsoid dimensions, which are then matched against the GRASS ellipse.table to find the corresponding GRASS ellipsoid name. I suppose we could add a list of GRASS ellipsoid names and PROJ equivalents somewhere as a way around this.

Datum names would also be a problem, but a very small one as PROJ only supports 10 named datums. One of these is not in GRASS and could be added, 5 of the others have exactly the same short name, and the other 4 differ only in the case of the letters (PROJ- upper case, GRASS- lower case).

So if an additional field was added to ellipse.table with a PROJ-compatible name, datum GGRS87 was added to GRASS, and PROJ datum names were converted to lower case before being imported into GRASS, it should be possible to write an OSR-independent function in lib/proj/convert.c to convert from PROJ to GRASS format (the OSR-dependent parts of lib/proj are already conditionalised on HAVE_OGR; there would be no problem in adding more functions outside of that).

Paul

10

On Thu, 21 Aug 2008, Paul Kelly wrote:

On Thu, 21 Aug 2008, Glynn Clements wrote:

Paul Kelly wrote:

AFAICT, the only "core" GDAL dependencies are the OSR stuff in
lib/proj and g.proj, and the vector library's support for v.external.
In which case, it really ought to be possible to get a usable version
of GRASS with no GDAL dependency in any of the libraries or core
modules (e.g. g.proj, g.region).

It used to be possible to compile g.proj without GDAL support, and use it
only for reporting on the projection information and verifying datum
information, but I removed all the conditional stuff for that when GDAL
became a non-optional dependency. Could be put back in I suppose.

That would be useful.

I'll get on to it but it won't be before the end of the weekend I'd say.

I've committed some changes to SVN trunk now to allow g.proj to compile with reduced functionality if HAVE_OGR is not defined.

Also, is there any fundamental obstacle to being able to set the
projection information using PROJ parameters without going via OSR?
Given that GRASS uses PROJ for the actual projection, it shouldn't
actually need anything other than PROJ parameters. g.setproj doesn't
use any lib/proj functions except GPJ_ask_datum_params() (which no
longer exists).

The main issue is the short ellipsoid names, which can be quite different. I believe the GRASS ellipse.table pre-dates the integration of PROJ into GRASS in the early 1990s. OSR understands the PROJ ellipsoid names, so at present it is used to get the underlying ellipsoid dimensions, which are then matched against the GRASS ellipse.table to find the corresponding GRASS ellipsoid name. I suppose we could add a list of GRASS ellipsoid names and PROJ equivalents somewhere as a way around this.

Datum names would also be a problem, but a very small one as PROJ only supports 10 named datums. One of these is not in GRASS and could be added, 5 of the others have exactly the same short name, and the other 4 differ only in the case of the letters (PROJ- upper case, GRASS- lower case).

So if an additional field was added to ellipse.table with a PROJ-compatible name, datum GGRS87 was added to GRASS, and PROJ datum names were converted to lower case before being imported into GRASS, it should be possible to write an OSR-independent function in lib/proj/convert.c to convert from PROJ to GRASS format (the OSR-dependent parts of lib/proj are already conditionalised on HAVE_OGR; there would be no problem in adding more functions outside of that).

Been thinking about this a bit more and a more fundamental problem than
the ellipsoid and datum names is the lack of validation of the PROJ.4
parameters that would be possible. All g.proj could do would be to pass
the parameters into the PROJ_INFO and PROJ_UNITS unchanged and any
problems would only show up when they were used for reprojection. I'm not
sure how good the error checking is there as it assumes that the PROJ_INFO
and PROJ_UNITS are valid. g.setproj always creates a valid file because it composes the parameters from a set list. g.proj currently uses OGR to parse and validate the PROJ parameters before converting them to a PROJ_INFO and PROJ_UNITS. So I feel that modifying g.proj to read a set of PROJ parameters and pass them directly into PROJ_INFO and PROJ_UNITS would lose a piece of functionality, i.e. validating the co-ordinate system. And it still would be complicated. At the minute I think it needs a bit more thought or discussion.

Paul

11

Paul Kelly wrote:

>> Also, is there any fundamental obstacle to being able to set the
>> projection information using PROJ parameters without going via OSR?

Been thinking about this a bit more and a more fundamental problem than
the ellipsoid and datum names is the lack of validation of the PROJ.4
parameters that would be possible. All g.proj could do would be to pass
the parameters into the PROJ_INFO and PROJ_UNITS unchanged and any
problems would only show up when they were used for reprojection. I'm not
sure how good the error checking is there as it assumes that the PROJ_INFO
and PROJ_UNITS are valid. g.setproj always creates a valid file because it
composes the parameters from a set list.

So, it looks like I need to cannibalise g.setproj so that there is
still a way to set projection parameters in the absence of OSR.

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