[raster] Two-raster spatial operations framework proposal

[bnordgren]

Note: this patch requires #1050 and #1053, and introduces a dependency from rt_pg to postgis-2.0.so. The patch applies against #7453.

I doubt this is usable yet, but there's enough here to start people talking (e.g., draw fire). This is a C based framework for two-raster operations returning a raster result. The two "major phases" of such a task are: precalculate the expected extent of the result; and iterate over the result raster to produce one value per cell. Both major phases have been encapsulated in their own function, and both use function callbacks as a mechanism to implement extension points. The framework takes differing alignments, cell sizes, and projections into account, such that the callbacks need only to worry about their fundamental task.

---

The grid extent function has one extension point defined: calculate extent of result. The remainder of the function retrieves the extents of the two input rasters and projects them into the result srid, creates a raster of the appropriate size, and sets various metadata fields. The four current implementations of this extension point calculate extents for "intersection", "difference", "union" and "symdifference".

---

The result-raster-iteration function has two extension points defined: determine whether a cell is included in the result; and provide a value for the cell.

There are four implementations of "determine whether a cell is included in the result", corresponding to the four extent calculations: "intersection", "difference", "union" and "symdifference".

There are two implementations of "provide a value for the cell": one generates a mask; the other copies image data from one of the input rasters (r1 if possible, r2 if not). The implementation which copies image data sets "not present" cells to nodata.

---

Provided one can access Postgresql's expression parser from C, it should be possible to implement a 2 raster MapAlgebra by implementing a "provide a value for the cell" extension. The interface may have to be augmented slightly to allow "user" parameters to be passed (such as the expression itself, as well as the bindings from symbol to band)

[bnordgren]

[bnordgren@ghosty-mosty postgis-raster]$ svn info
Path: .
URL: http://svn.osgeo.org/postgis/trunk
Repository Root: http://svn.osgeo.org/postgis
Repository UUID: b70326c6-7e19-0410-871a-916f4a2858ee
Revision: 7453
Node Kind: directory
Schedule: normal
Last Changed Author: pramsey
Last Changed Rev: 7453
Last Changed Date: 2011-06-23 15:36:54 +0000 (Thu, 23 Jun 2011)

[bnordgren@ghosty-mosty postgis-raster]$ svn status
M       raster/rt_pg/Makefile.in
M       raster/rt_pg/rtpostgis.sql.in.c
M       raster/rt_pg/rt_pg.c
M       raster/rt_core/rt_api.c
M       raster/rt_core/rt_api.h
M       configure.ac

[bnordgren]

Linking to the shared postgis library in the build directory is causing problems at runtime (after install). I could use some help on this one:

psql:rtpostgis.sql:39: NOTICE:  type "raster" is not yet defined
DETAIL:  Creating a shell type definition.
psql:rtpostgis.sql:39: ERROR:  
  could not load library "/usr/lib/postgresql/rtpostgis-2.0.so": 
   /home/bnordgren/build/local/postgis-svn/src/trunk-build/postgis/postgis-2.0.so:
     cannot open shared object file: Permission denied

[mcayland]

I think that the way in which you're doing this is wrong and likely to cause strange runtime errors depending upon the linkers in different platforms.

I would suggest that you submit as 2 separate patches: the first to abstract and move the PROJ.4 cache API into liblwgeom, and then a second implementing your new functionality.

[bnordgren]

Clearly. :)

I can certainly submit two patches. At the abstract level, does this mean that all common code (e.g., required by both ./postgis and ./raster) needs to be factored out to liblwgeom and there can be no dependencies on postgis-2.0.so by other modules?

I was just a little leery of being too heavy handed moving things around for fear of irritating people.

Will break this out over the weekend.

[mcayland]

Excellent. For extra bonus points you could rework the caching API so that it can be reused in multiple places; for example I believe both prepared geometries and PROJ.4 use a per-query cache and so it would be nice to unify this code and make it easily available to other functions.

[bnordgren]

Ok, this is what I get for opening my mouth.

It appears to be more complicated than I thought at first blush (or so says the compiler). The bit I need to depend on is the proj4 cache api. This part connects to the database to get the proj4 string from the spatial_ref_sys table. Whoops. It also depends rather heavily on the postgresql memory context. Whoops. Whoops.

Moving this code to liblwgeom would create a dependency from liblwgeom to the core postgresql libraries, and would probably make us compile liblwgeom as a postgresql module, if I understand this pgxs thing right.... I need to reconfirm with you that this is actually what you want to do before I up and do it.

[mcayland]

:) Well let me clarify how I would see things working.

Firstly: liblwgeom is intended to be standalone in that it should be possible to make use of it without PostgreSQL present at all. So the current route of trying to link against postgis.so is incorrect so please don't spend any more time on it.

I'd envisage something along the following lines:

1) Create a new file liblwgeom/lwgeom_proj.c

2) Add some basic APIs which take a proj4text and return a projPJ handle

3) Move the main reprojection function from postgis/lwgeom_transform.c into liblwgeom/lwgeom_proj.c taking a projPJ handle as a primary parameter

4) Rework postgis/lwgeom_transform.c so it calls the underlying liblwgeom routines instead

5) Rework postgis/lwgeom_transform.c cache code so it uses the underlying liblwgeom routines... and while you're there, also refactor it into a simple API so that adding caching to any other functions within PostGIS (as in linked with PostgreSQL) becomes much easier (and indeed, I believe you can test this by making the prepared geometry cache use your new code too).

Does that make more sense now?

Mark.

[bnordgren]

Replying to mcayland:

Firstly: liblwgeom is intended to be standalone in that it should be possible to make use of it without PostgreSQL present at all. So the current route of trying to link against postgis.so is incorrect so please don't spend any more time on it.

I think I was unclear. Here's the situation:

Liblwgeom depends on external libraries, such as geos. [OK]

Raster depends on liblwgeom [OK]

Raster also depends on postgis now [PROBLEM].

It's the link between rtpostgis-2.0.so and postgis-2.0.so (with raster depending on postgis) which is problematic. In terms of your outline, this is the part which requires that problematic link:

5) Rework postgis/lwgeom_transform.c cache code so it uses the underlying liblwgeom routines... and while you're there, also refactor it into a simple API so that adding caching to any other functions within PostGIS (as in linked with PostgreSQL) becomes much easier (and indeed, I believe you can test this by making the prepared geometry cache use your new code too).

Do you see a way to avoid making this link?

[mcayland]

Replying to bnordgren:

I think I was unclear. Here's the situation:

Liblwgeom depends on external libraries, such as geos. [OK]

Raster depends on liblwgeom [OK]

Raster also depends on postgis now [PROBLEM].

It's the link between rtpostgis-2.0.so and postgis-2.0.so (with raster depending on postgis) which is problematic. In terms of your outline, this is the part which requires that problematic link:

5) Rework postgis/lwgeom_transform.c cache code so it uses the underlying liblwgeom routines... and while you're there, also refactor it into a simple API so that adding caching to any other functions within PostGIS (as in linked with PostgreSQL) becomes much easier (and indeed, I believe you can test this by making the prepared geometry cache use your new code too).

Do you see a way to avoid making this link?

I think so? I was suggesting refactoring the reprojection code and moving its core into liblwgeom. Since raster already depends upon liblwgeom, then the problem is already solved and all other users of liblwgeom can also benefit from the revised functionality.

[bnordgren]

Replying to mcayland:

5) Rework postgis/lwgeom_transform.c cache code so it uses the underlying liblwgeom routines... and while you're there, also refactor it into a simple API so that adding caching to any other functions within PostGIS (as in linked with PostgreSQL) becomes much easier (and indeed, I believe you can test this by making the prepared geometry cache use your new code too).

Do you see a way to avoid making this link?

I think so? I was suggesting refactoring the reprojection code and moving its core into liblwgeom. Since raster already depends upon liblwgeom, then the problem is already solved and all other users of liblwgeom can also benefit from the revised functionality.

There may be a couple of problems here. First, the cache seems to make use of "MemoryContext", which I believe is a postgresql structure. Presumably, this is to ensure that things which are allocated are free-ed at the end of the call. However, they seem to be storing at least a little state information in there too.

Second, raster needs access to more than just the "core code". It also needs access to the code that retrieves the proj4 string from the spatial_ref_sys table. If we can't find a way to link raster to postgis, we'll have duplicate code in both libraries.

[pracine]

Bborie wrote something to get the projection info in RASTER_transform(). Did you have a look? Isn't that doing what you need?

[bnordgren]

Replying to pracine:

Bborie wrote something to get the projection info in RASTER_transform(). Did you have a look? Isn't that doing what you need?

The projection stuff is for LWGEOMs (essentially, the "extent"s of the input rasters.)

I'm actually not projecting rasters. When I need to access a raster from a different projection, I'm using GDAL to calculate the indices in one image given the indices in another. Check it out:

​http://www.gdal.org/gdal__alg_8h.html#7671696d085085a0bfba3c3df9ffcc0a

[Bborie Park]

You may want to take a look at rt_raster_gdal_warp as that provides most of the functionality provided by the gdalwarp utility.

[bnordgren]

Blocked by #1133. Much of the discussion of how specifically to divide functionality between libraries applies to #1053 instead of this ticket (although #1053 is a dependency of this ticket).

[bnordgren]

Most recent patch reflects "option 3" on #1133. It applies against r7683 which has already had patches from #1050 and #1053 applied.

This patch:

• Moves the PG_MODULE_MAGIC and lwgeom_init_allocators to libcommon. (The latter uses the main postgis memory management functions.)
• Leverages the libcommon and liblwgeom static libraries.
• Requires the GEOS operations on the LWGEOM type (moved to liblwgeom in #1050)
• Requires the transformation API cache used by postgis (moved to libcommon in #1053).

All cunit tests pass except for loader/Latin1, which fails on all instances of r7683 on my platform. Regression tests pass. I did not see it run additional regression tests for raster (only for postgis).

[bnordgren]

Updated patch applies against r7720 with #1053 applied. It is no longer necessary to apply #1050 separately.

[bnordgren]

Updated patch against r7728

[robe]

Milestone renamed