[raster] ST_Transform

[dustymugs]

This function would provide rasters the ability to do on-the-fly reprojections. It may be possible that underlying code could be shared with the planned ST_Warp and ST_Resample functions but for now, the code will focus on supporting ST_Transform.

A significant limitation with the C warp API in GDAL is that it is a subset of the C++ API. Based upon my review of the C warp functions available and the code for  gdalwarpsimple.c, it looks like we are severely constrained by the function GDALSimpleImageWarp. According to the docs, GDALSimpleImageWarp's "algorithm is called simple because it lacks many options such as resampling kernels (other than nearest neighbour), support for data types other than 8bit, and the ability to warp images without holding the entire source and destination image in memory."

So, my question is, do we go C++, put ST_Transform on the backburner while the C API is improved or do what we can with the current C API?

[pracine]

Could we ask on the GDAL list why only a limited version is available in the C API and what would be the effort of making the full version available? We could takes a wiser decision.

[dustymugs]

Yes. We'll probably have to. I'll post a message later tonight. Any thoughts on what to do in the interim? ST_Transform is a critical enough function that I wouldn't mind making minimal use of the C++ API though a quick audit of the PostGIS codebase shows a strict adherence to only using C.

[pracine]

So what would you suggest? A C++ shared library exposing some C functions (ideally similar to the ones GDAL should expose)? Do we actually need a shared library? Not just a C file exposing C++ methods? I'm not really familiar with C++ binding from C.

[dustymugs]

I don't think we'd want to maintain an additional shared library interfacing between the C++ and C. I was thinking in terms of directly contributing patches to GDAL so that we can continue the current paradigm of sticking with C and being dependent on GDAL.

I plan on doing some work the rest of this week on GDALPolygonize to add support for floating points and rounding, so that may give me some ideas of what needs to be done to extend GDALSimpleImageWarp. Like you, I've never worked on making C bindings from a C++ API. But if we wait for someone to develop it in GDAL, we may be waiting forever at the expense of PostGIS Raster.

[pracine]

:-) I'm very happy with that! Vive l'open source!

[dustymugs]

For PostGIS 2.0, a basic ST_Transform function should be provided. These functions use GDAL's Warp API, specifically the C function GDALAutoCreateWarpedVRT(). This function allows us to quickly provide an ST_Transform function for rasters at the expense of not doing multiple actions at the same time as the reprojection (e.g. changing the pixel scale or the skew). It is expected that the C function GDALAutoCreateWarptedVRT() will be replaced with a more generic function matching that used in the GDAL utility gdalwarp.

1. ST_Transform(rast raster, srid integer[, algorithm text DEFAULT 'NearestNeighbor?']) -> raster

raster: the raster to reproject to a new spatial reference system. if raster's SRID is -1, the function returns NULL

srid: the id of the spatial reference system to reproject the raster

algorithm: the resampling algorithm to use. This argument is case-insensitive. The default algorithm is "NearestNeighbor?". Possible algorithms are

NearestNeighbour

Bilinear

Cubic

CubicSpline

Lanczos

The docs for GDAL function GDALAutoCreateWarpedVRT() does not indicate if the Lanczos algorithm is supported but is provided above in the list of possible algorithms for the eventual replacement of the GDALAutoCreateWarpedVRT() function with a more generic GDAL warp function

ST_Transform(raster, 3310)

ST_Transform(raster, 4326, 'Cubic')

ST_Transform(raster, 3309, 'Bilinear')

[dustymugs]

I need to amend the number of function parameters for ST_Transform. The correct version is:

1. ST_Transform(rast raster, srid integer[, algorithm text DEFAULT 'NearestNeighbor'[, maxerr double precision DEFAULT 0.125]]) -> raster

maxerr: max error allowed in approximating the transformation (0.0 for exact, no approximations). GDAL defaults to 0.125.

[darkblueb]

Replying to dustymugs:

I need to amend the number of function parameters for ST_Transform.

Hey - sorry to show up late, but, please confirm that this has no effect on the existing st_transform() call used by soooo much code out there, for geometry

[dustymugs]

I can confirm that this won't affect anything related to the ST_Transform function for geometries. This is purely an analog to the geometry version for rasters.

[dustymugs]

Added in r7338.

[pracine]

Great!

What do you mean by "at the expense of not doing multiple actions at the same time as the reprojection (e.g. changing the pixel scale or the skew)"?

Do we need a special build of GDAL?

[robe]

Just a note -- I'm sure I just have to tweak my config but crashes on testapi with libproj-0.dll not found (couldn't load library). Why it's looking of libproj-0.dll instead of my libproj.dll I don't know.

[dustymugs]

Replying to pracine:

Great! What do you mean by "at the expense of not doing multiple actions at the same time as the reprojection (e.g. changing the pixel scale or the skew)"? Do we need a special build of GDAL?

You don't need a special build. The standard build is just fine. When I use the utility gdalwarp, I like to specify the pixel size of the output raster as the pixel size can and probably will change with the reprojection. The current underlying GDAL function doesn't provide the opportunity to specify the GeoTransform? matrix as it autogenerates it. The code for gdalwarp shows how to do it but I haven't had the time to walk my way through it.

[dustymugs]

Replying to robe:

Just a note -- I'm sure I just have to tweak my config but crashes on testapi with libproj-0.dll not found (couldn't load library). Why it's looking of libproj-0.dll instead of my libproj.dll I don't know.

I think you will need to tweak your config. Is it possible that your GDAL library wasn't built against your Proj.4 library?

[robe]

That's what I was thinking since I didn't add any extra parameters. I'll try to rebuild it pointing it at my proj config and see if that fixes the issue. If so I guess we'll need to note that in the docs that you have to build and may need to specify your proj location. The other thing to keep in mind at least for windows is that we have a somewhat hacked proj config that looks for proj in the share/contrib/postgis/proj folder. As I recall I don't think this is specifically for windows though I think Mark had done it with that in mind.

Ideally we want it using the same proj epsg/shift files and proj library as the rest of PostGIS, otherwise I can see all sorts of issues with differences in proj coming up. I already have such issues with my Linux box and windows box using a different version of proj library (or at least I think that is part of my headache with non-related to this things).

I didn't check does it use the spatial_ref_sys table or rely on the epsg files in proj folder? As I recall I think PostGIS just uses the datum shift files on not the other files.

Anyrate this is great work. Can't wait to try it out.

[dustymugs]

This ST_Transform function uses the values in the spatial_ref_sys table as I thought that'd be the most consistent implementation across platforms. I didn't know about the epsg files in the proj folder though I'll probably never use it anyways. GDAL covers up the gory details about what is needed for proj.4. All I do is pass along the srtext for the input raster and the output SRID.

[robe]

I think we still need to worry about the datum shift files though since even PostGIS relies on those. Hopefully Mark C. will chime in on this. He was the one that put in logic in PostGIS for it to look in specific place and then fall back on default setting (or is it vice versa), can't remember.

I think the code is on this line: http://trac.osgeo.org/postgis/browser/trunk/postgis/lwgeom_transform.c#L602

[dustymugs]

Good point. I have no idea has GDAL does it. I'll have to test by reprojecting some raster into NAD83...

[robe]

Pick a raster that's NAD27. NAD27 to NAD83 I think is when you run into issues.

[dustymugs]

Will do later today. I'm currently working on adding coverage table support to ST_Histogram.

[pracine]

I also have the "Unable to load PROJ.4 library (libproj-0.dll)" message. So what's the trick?

[robe]

I have no idea -- I didn't see a configuration option in GDAL to allow me to specify the proj library (I recall that when I compiled proj 4.7 it did name the file libproj-0.dll).

I went ahead and tested the examples in my PostgreSQL 9.1 install. It doesn't give an error but returns the wrong answer. Maybe new gdal expects proj 4.7 or above though it should error out and its not. I'm running proj 4.6.1.

So via psql -

SELECT
	ST_SRID(rast),
	ST_Width(rast),
	ST_Height(rast),
	ST_NumBands(rast),
	ST_ScaleX(rast),
	ST_ScaleY(rast),
	ST_Count(rast),
	ST_Mean(rast)
FROM (SELECT ST_Transform(
	ST_SetValue(
		ST_SetValue(
			ST_SetValue(
				ST_AddBand(
					ST_MakeEmptyRaster(10, 10, -500000, 600000, 1000, 1000, 0, 0, 1002163)
					, 1, '64BF', 0, 0
				)
				, 1, 1, 1, -10
			)
			, 1, 5, 4, 0
		)
		, 1, 5, 5, 3.14159
	)
	, 1003310) AS rast
) foo;

Give me this -- fair enough right:

 st_srid | st_width | st_height | st_numbands | st_scalex | st_scaley | st_count
 |  st_mean
---------+----------+-----------+-------------+-----------+-----------+----------+-----
 1003310 |       10 |        10 |           1 |      1000 |     -1000 |        2 | -3.429205

Note in BBorie's regress, he gets a width of 12 , scale x and scale y changed and I'm getting 10 back like its not doing anything but not erroring out either.

[robe]

actually just noticed that BBorie's examples would involve using the datum shift files which maybe it's not finding in my install. I'll try with ones that don't require datum shifts to see if it produces something more realistic.

Bborie, maybe you can add some that don't require datum shift like say NAD83 california feet to NAD83 long lat (4269).

[dustymugs]

I've added a test that don't require datum shifts from 3310 to 4269 in r7343.

I am using proj-4.7.0 on all of my dev/test machines.

[dustymugs]

It is possible (and likely if our experience with ST_AsGDALRaster w/ various GDAL version is repeated) that the output you're getting is correct. It isn't ideal and may force me to decide on what would be an appropriate test, if one consistent test is even possible.

[robe]

I'm using a fairly recent GDAL trunk -- 1.9 SVN r22508 (so that I could get the changes that you and Jorge made to GDAL)

Tried one of your new ones. Sadly that doesn't seem to work either -- so it seems it's not transforming on mine but silently failing and returning the same thing.

SELECT
	ST_SRID(rast),
	ST_Width(rast),
	ST_Height(rast),
	ST_NumBands(rast),
	round(ST_ScaleX(rast)::numeric, 3),
	round(ST_ScaleY(rast)::numeric, 3),
	ST_Count(rast),
	round(ST_Mean(rast)::numeric, 3)
FROM (SELECT ST_Transform(
	ST_SetValue(
		ST_SetValue(
			ST_SetValue(
				ST_AddBand(
					ST_MakeEmptyRaster(10, 10, -500000, 600000, 1000, 1000, 0, 0, 1003310)
					, 1, '64BF', 0, 0
				)
				, 1, 1, 1, -10
			)
			, 1, 5, 4, 0
		)
		, 1, 5, 5, 3.14159
	)
	, 1004269) AS rast
) foo;

Gives pretty much what I started with:

1004269	10	10	1	1000.000	-1000.000	2	-3.429

I'll try compiling proj 4.7 and see if it behaves differently or it still needs those files for some reason. Perhaps I'll set the proj data path on my box to see if it changes the behavior first. I have been putting off compiling proj 4.7 just because it's cumbersome to compile under mingw windows -- always have to hack it (come to think of it only been succesful compiling it under my mingw-64 install though haven't tried for a while on my regular mingw.

[dustymugs]

As you do that, I'll recompile my stack for proj 4.6.1 and see how the tests go in linux.

[dustymugs]

Running proj 4.6.1, gdal 1.9 r22520 and postgis r7343, all regression tests pass on a 32-bit linux box. Will test 64-bit linux box next.

[dustymugs]

Just finished testing on a 64-bit linux box running same stack (proj 4.6.1, gdal 1.9 r22520 and postgis r7343) as the 32-bit version. All regression tests pass.

I don't know if I want to go into the Windows world...

[robe]

Nah you don't want to go into the windows world -- no one does -- except for us though :)

It's a very dark world with software license monsters lurking in every corner. strk will attest to that. Seriously though it could be I'm compiling gdal wrong.

Is there a switch to specify the proj to link to like we have for postgis configure? I thought maybe it was relying on the GDAL_DATA and PROJ_LIB paths and so forth at runtime, but I added those to my environment variables and still doesn't seem to make a difference. It's possible I have a typo somewhere though.

[dustymugs]

I don't believe there is a switch to specify the proj. I think I'm going to have to bite the bullet and install Windows (probably Windows XP or 7, suggestions?) so that I can see the problems you're having.

[robe]

Bborie,

Either one should work fine. I have an old windows xp and a new windows 7 and have been able to compile fine and run on both. So whatever is most accessbile. Windows xp is a bit less finicky with security stuff so often eaiser to get going without changing permissions, but that's about it.

FWIW I tried adding my proj to my search path and recompiling GDAL. I was also able to compile GDAL trunk with libtool on (so guess they must have fixed that issue). Neither fixed the issue unfortunately. I also turned on filemon to see if while I was running the transform it was looking for a path, but didn't see it looking for anything aside form the files that postgres usually loads up.

Though I didn't recompile postgis after I recompiled gdal, so I'll try that next.

Pierre -- Are you having the same issue and are you running on Windows or your Linux box. Before Bbories goes on a wild goose chase and has to suffer the torture of trying to get things working on windows, perhaps we should confirm that it's only a windows problem and no other Linux folks are having the issue?

[pracine]

I'm still and always on Windows. But Windows 7 now. I searched for libproj*.dll on my system and it doesn't exist so I guess PostGIS link statically with it when PostGIS raster is linking dynamically with it. Make sence?

Bborie I don't recommend trying setting you up on Windows. You'll lose days...

[mcayland]

Hmmm so do you have a static or a dynamic PROJ.4 link on Windows? The trick to producing the DLL is this line on http://trac.osgeo.org/postgis/wiki/DevWikiWinMingWSys_14_15:

gcc -shared -o libproj.dll -Wl,--out-implib=libproj.dll.a -Wl,--export-all-symbols -Wl,--enable-auto-import -Wl,--whole-archive libproj.a -Wl,--no-whole-archive /c/mingw/lib/libmingw32.a

Perhaps you need to temporarily rename libproj.a after you've done this in order to prevent raster finding the static version?

In terms of the PROJ.4 search paths for grid files, Regina mentioned this:  ttp://trac.osgeo.org/postgis/browser/trunk/postgis/lwgeom_transform.c#L602. All this does is ADD this to the search path, so PROJ.4 will still search in other locations. I believe that these are currently the value of the PROJ4LIB environment variable (if it exists), then the specified search path, and then finally a fixed location generated at compile time (typically C:\PROJ\LIB).

Hope this helps you guys :)

[dustymugs]

You two are really trying to dissuade me from dealing with the ordeal of a Windows build environment ;-). But, I've taken the plunge and am building a Window XP VM. From the looks of it, is everything compiled through MinGW? Any reason we don't do Visual Studio or something of that sort either (since GDAL officially doesn't support MinGW)?

I always like have some approximation of what end-users are using so that when these kinds of problems come up, I'm not flailing my hands wildly...

[pracine]

GDAL do not support MinGW officially but we, after some struggle, succeeded to build it last winter. So I'm still using libgdal.dll 1.7.1

PostGIS do not support VC yet so you have no choice than tweaking a build.

On Windows you are always condemned to tweak things and it can take days...

I'm still struggling with my build. My understanding is that GDAL is looking for the wrong a version of libproj. Even if I build a dll like Mark suggested, it does not use it. No more success if I rename it libproj-0.dll.

GDAL 1.8 do not compile either. So I try to make my 1.7.1 to work. I really don't want to throw everything to the garbage and start from scratch. Really a nice day...

[robe]

see #1015 for more detailed example on Linux

[dustymugs]

Can someone with the problem try r7364? I've explicitly set the option that is causing the error message posted by strk in #1015.

[pracine]

Tested r7365. Same problem:

ERROR 6: Unable to load PROJ.4 library (libproj-0.dll), creation of OGRCoordinateTransformation failed. GDALDataset Check failed on line 1490 make[2]: *** [check] Error 1 make[2]: Leaving directory `/c/thesrc/postgis/postgis-trunk/raster/test/core' make[1]: *** [core-check] Error 2 make[1]: Leaving directory `/c/thesrc/postgis/postgis-trunk/raster/test' make: *** [core-check] Error 2

[dustymugs]

Interesting. That detailed error message tells me that the problem I fixed is something different. strk, can you try again with r7364 or later?

Pierre, what I can tell you based on the "creation of OGRCoordinateTransformation failed" part is that you'll need to set the environment variable PROJSO. See:

 http://trac.osgeo.org/gdal/wiki/ConfigOptions#PROJSO

[pracine]

Setting PROJSO works... Thanks BBorie.

Now the rt_transform.sql regress fail:

--- rt_transform_expected	2011-06-10 14:45:17 -0400
+++ /tmp/pgis_reg_5856/test_45_out	2011-06-10 15:19:16 -0400
@@ -4,9 +4,9 @@
 1003309|12|12|1|995.277|-995.277|2|-3.429
 1004269|12|9|1|0.011|-0.011|1|-10.000
 |||||||
-1003310|12|12|1|995.262|-995.262|7|-1.427
+1003310|12|12|1|995.262|-995.262|5|0.513
 |||||||
-1003310|12|12|1|995.262|-995.262|20|-0.515
+1003310|12|12|1|995.262|-995.262|5|0.513
 1003310|10|10|1|1000.006|-1000.006|2|-3.429
 1003309|10|10|1|999.994|-999.994|2|-3.429
 BEGIN

I have GDAL 1.7.1

[robe]

okay will give this a try too.

[robe]

We do have to think about having it look for this somehow though like how postgis looks for it in a certain place. Bborie, is that possible or is it too deep in GDAL code?

when installing PostGIS via StackBuilder?, I think it would be difficult to configure the script to put this in without possibly risking screwing up other peoples existing installs GDAL related stuff.

[dustymugs]

This is quite deep in GDAL code. The GDAL ticket #1338 discusses this.

 http://trac.osgeo.org/gdal/ticket/1338

[robe]

Well I'm not getting the error anymore abot libproj-0.dll, but I get check failed line 1490.

Check failed on line 1490
make[2]: *** [check] Error 1
make[2]: Leaving directory `/c/projects/PostGIS/trunk/raster/test/core'
make[1]: *** [core-check] Error 2
make[1]: Leaving directory `/c/projects/PostGIS/trunk/raster/test'
make: *** [core-check] Error 2

and still get the wrong answers in postgresql. Pierre -- are you answers right in PostgreSQL? If you just run the regress tests in an installed raster db?

As far as the proj thing given that it seems to default to libproj-0.dll or some such thing, which I recall is the name I get when I compile proj 4.7.1, it might not be that bad.

I was thinking of releasing PostGIS 2.0 with proj 4.7 anyway, but wanted to verify it didn't break things -- which as I said I'm having some issues with my linux box installed with prj 4.7 not being able to reproject things which is part reason I can't use my Linux box for production work. But that's neither here nor their, but I guess testing on 4.7 on windows will rule out issue only on linux with proj differences.

[robe]

Bborie,

False alarm -- the tests pass if I add my project 4.6.1 to my path settings.

Though I think I sttil have issue geting right answers under PostgreSQL. That might be something different.

[robe]

Bborie,

It works under PostgreSQL now. I thought I had my global environment PROJSO set, but I had a typo. So fixing that worked. Though I could have sworn the first time they all registered as NULL -- all the fields. But can't recreate that again. So I guess you don't need to venture into windows land.

strk -- can you confirm things work for you with PROJSO environment variable then I guess we can close this out.

[dustymugs]

Pierre,

It's quite strange as to what is happening in the regression. When I test on 1.7.3, the same two tests fail but with different values. Both tests are for different resampling algorithms (Bilinear and Cubic).

Regina,

What version of GDAL are you using?

[strk]

It's crowded here. I'll be in bug #1015.

[robe]

not absolutely sure -- gdal trunk no older than a day or 2 ago.

[robe]

BBorie,

Sorry I think I'm still having problems all the noise from my 9.1 build that always fails because of that issue with ST_MapAlgebra got me confused and also I think some of the tests past because I had put in some extra gdal environment variables like gdal-data, proj_lib and when I thought things were working I took them out.

I did see 12s in there so I suspect one of those other env variables helped before. I'm haunting that down.

[dustymugs]

Pierre and Regina,

Have either of you been able to get the regression tests to run successfully? I did make some changes to the regression tests due to the different answers possible if using different versions of GDAL.

If neither of you have problems, can one of you close this ticket? Thanks!

[robe]

Bborie,

My regress test on 8.4 passed fine (the transform check), minus the 2 non transform failures you probably are also getting because of the GSERIALIZED changes. I didn't close this out yet though because I wasn't sure if it was the path settings I had put in or the change in your regress tests. So wanted to pull my path settings out before I close.

Also that ST_MapAlgebra issue is making it difficult for me to test on my 9.1 and 9.0 as most tests fail on those.

[dustymugs]

I just committed a patch r7396 which resolves the regression failures. The regression tests were trying to use SRID values that exceeded the allowed maximum of GSERIALIZED.

[robe]

All my raster tests seem to be passing now. Though I guess the real test will be when I combine ST_Transform with ST_AsJPEG etc.

[pracine]

What is the meaning of the 'maxerr' parameter?

[dustymugs]

The maxerr parameter is the "error threshold for transformation approximation (in pixel units - defaults to 0.125).". It is the same as the "-et" parameter of gdalwarp.

 http://www.gdal.org/gdalwarp.html

[pracine]

...and for my grand mommy?

[dustymugs]

Replying to pracine:

...and for my grand mommy?

???? I don't get it.

If you don't want any approximation, you can set "maxerr" to zero (0.). I set the default to the same as gdalwarp.