transform() grid-shift 2nd chance logic defective

[william....@…]

What steps will reproduce the problem?

select x(transform(pointfromtext('POINT (173 52)',4269),4267)), y(transform(pointfromtext('POINT (173 52)',4269),4267)); WARNING: transform: -38 (failed to load NAD27-83 correction file) ERROR: transform: couldn't project point: -38 (failed to load NAD27-83 correction file) <SELECT aborted; no output rows>

What is the output as corrected?

select x(transform(pointfromtext('POINT (173 52)',4269),4267)), y(transform(pointfromtext('POINT (173 52)',4269),4267)); WARNING: transform: -38 (failed to load NAD27-83 correction file) WARNING: transform: -38 (failed to load NAD27-83 correction file)

x | y

---

173 | 51.9999999990852

(1 row)

Please use labels and text to provide additional information.

The PostGIS transform() function contains logic that intends to substitute a '2nd chance' non-shifted transformation for an unavailable grid-shift transformation (e.g., NAD83 to NAD27 transformation of a point outside of the area covered by the NAD27 grid-shift table, or in the absence of that table). The logic handling this situation is located in the transform_point() function of lwgeom_transform.c but seems outdated with respect to how recent versions of proj.4 handle the situation.

This discord causes a 'SELECT … transform(geom,srid) … FROM mytab' to abort and be rolled back where even a single row of table mytab contains a geometry geom that lies outside the bounds of an applicable grid-shift table. Running on Postgres 8.3.5, PostGIS 1.3.3, and proj.4 4.6.1, such a SELECT aborts with the twin messages shown in the problem reproduction above.

Attached is a patch that seems to be a reasonable solution for this problem.

Although this patch works for me, please scrutinize it carefully. This is my first look at PostGIS and proj.4 source code. Worse, I code in C very rarely and am especially error-prone when using C pointers and indirection.

I have not been able to fabricate data that exercises every branch of this

patch, so it may still contain errors.

Thanks so much for PostGIS! It's a tremendous resource that contributes materially to the public good.

[mcayland]

Hi William,

Thanks for submitting the patch with your bug report above. My main question is that you state 'The logic handling this situation is located in the transform_point() function of lwgeom_transform.c but seems outdated with respect to how recent versions of proj.4 handle the situation' but the new method of handling is not clear - can you point towards some documentation?

I'm also a little surprised by the existing logic to try again without the grid shift information - surely we should simple return NULL in the case where the transformation is invalid?

Also another minor note: if you submit a patch, could you submit in either unified (-u) or context (-c) format? It just makes it easier to review on sight without application.

ATB,

Mark.

[william....@…]

Mark,

Thanks for the prompt reply. Didn't really expect anyone to look at this until 2009.

Answers to your questions …

1. Documentation re. new method of handling issue? First, I noticed that proj's

'cs2cs' utility does NOT throw an error (or even a warning) in this situation — rather, it simply performs the transform (which is no transform in this case) without the grid shift. I then read cs2cs's relevant code and pretty much plagiarized it for my patch.

Three other supporting elements:

(a) lwgeom_transform.c's '2nd chance' remark in function transform_point().

(b) lwgeom_transform.c's function pj_transform_nodatum(), which does not function as claimed its leading comments.

© Very helpful remarks for rev. 1.11 in the revision history contained in the leading comments of proj's pj_transform.c.

2. Why not return null when the transformation is invalid? In at least this case,

NAD83 → NAD27 transformation, the most useful behavior is to return an unshifted transformation. I have a large table of broadcast radio station locations — most in CONUS but some in Alaska and Puerto Rico; aborting processing of the whole table because the transformation of one or two rows will be off by a few hundred meters is very inconvenient, especially because the offending geometry is not identified in the error message (and only the first one is found before abort).

My main argument: if proj.4 thinks it's useful to return a non-shifted transform, PostGIS should faithfully propagate that behavior. With my patch, the user is warned of any 'transient' error (proj's terminology in pj_transform.c remarks). In this particular case, the warning is misleading ('failed to load NAD27-83 correction file' when the real problem is out-of-bounds input), but at least they're warned once for each offending geometry while still getting a useful result table. If you know how to make that warning actually identify the offending point, that would be a great improvement. Didn't want to press my luck given my very limited C skills.

3. Wrong diff format. I thought I'd seen instructions for the preferred form of

'diff,' but it looks like in updating PostGIS support Web pages to point to code.google.com, those instructions disappeared … or at least became impossible for me to find. I'm now attaching a 'diff -u'.

Thanks again for the quick response. Merry Christmas!

[mcayland]

Hi William,

I understand that in your case of a NAD83 → NAD27 transformation, ignoring the datum is an acceptable error to you, but I can't see this being the case with all applications and all datum shifts. We see people being caught out a lot because they've installed the grid shift files on one server (and not another) and wonder why the results beween 2 servers are different.

From a strictness point of view, I would argue that any transformation requested by the user that cannot be achieved using the current installation should throw an ERROR rather than trying to second-guess the user. At least then it would be very clear that they need to re-install PROJ.4 rather than wondering why they are getting erroneus results.

Frank, what are your thoughts on this?

ATB,

Mark.

[william....@…]

Mark,

Yep, I've had exactly the problem you cite: the Ubuntu and Debian PostGIS packages I was using lacked the grid-shift files and I was scratching my head wondering why my NAD83→NAD27 transformation did nothing. Had to build my own PostGIS to fix it — and that was long before I patched PostGIS for this issue. However, I believe my patch provides adequate safeguards for this problem. I believe the other reasons I give for the behavior I advocate are of equal or greater importance. To briefly summarize them:

• My patch throws a warning, so the user knows something is wrong unless they have

actively suppressed warnings. (Improving that warning by identifying the offending geometry would be an important improvement to my patch. I'll try to do it if the team doesn't have the resources to do so, though your C skills are surely far stronger than mine.)

• In all other cases (AFAIK), PostGIS simply 'wraps' proj.4 behavior and a

knowledgeable user expects to see that behavior propagated without change. PostGIS does not abort transformations outside the defined grid-shift region; its cs2cs utility glides past this problem without even a warning.

• I'm just a newbie to cartographic grid shifts, but from reading the proj.4

documentation, I don't think that NAD27 is unusual in being defined over a limited region. Anyone who's using grid shifts should know the limitations of what they're doing. Thus, if they present geometries for transformation outside the defined region, they presumably wish to obtain the transformation in all cases, even if it cannot be grid-shifted. I believe this behavior manifests the 'principle of least surprise.'

Whatever y'all decide, thanks so much for what you're doing!

— Bill

[william....@…]

Correction: 2nd bullet above, second sentence, first word, 'PosgGIS' should read 'proj.4'. Sorry.

[william....@…]

Mark, my previous comment omitted (a) the real root cause of this problem and (b) two possible fixes that might leave us both reasonably satisfied.

Real root cause, highlighted by your last remark: proj.4's error return code currently makes no distinction between (1) unavailable grid shift files and (2) points that fall outside the boundaries of an available grid shift file.

Best possible fix: PostGIS would throw an error in case 1 and throw a warning in case 2 — but there's presently no way for PostGIS to distinguish between the two cases. Any chance you could get the proj.4 team to add a distinct error code and message for case 2?

Second-choice fix: return NULL when proj.4 reports a grid-shift error (in either error case above). The user who wishes to obtain a non-grid-shifted transform as a fallback could then do so with coalesce(). On the other hand, a user who unexpectedly has rows with geometries that fall outside the area for which the grid shift is defined could easily identify all such rows in a single SELECT, which is impossible if you throw an error on the first offense. Finally, users who are missing a necessary grid-shift file at least won't get SLIGHTLY different answers from systems that have the necessary files; this should make it easier for them to find their real problem.

Throwing a generic error like PostGIS now does makes it very difficult for a user whose table contains a subset of rows with geometries outside a grid-shift region to identify those offending rows. The user must discover the boundaries of the grid shift file (how?), then craft a query to find geometries that fall outside tha area.

One last thought about improving the warning in my preferred fix: my current patch throws a warning for every offending POINT; this increases error-handling load by several orders of magnitude for complex geometries. Better would be to throw one warning for every GEOMETRY with an offending point, identifying the offendor by the first line of the geometry's summary() plus the coordinates of the first point in the geometry and the first offending point (if different). I haven't dug into the PostGIS transform code enough to see how difficult this might be.

Sorry to belabor this, but the current generic error behavior really cost me a lot of time, and I'd like to save others the inconvenience.

[fwarmerdam]

It is my opinion that a more appropriate location to adjust behavior is the +nadgrids specification for the coordinate system.

The default +nadgrids setting for NAD27 is (with +datum=NAD27):

+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat

The @ prefix means no error is reported if the files are not present, but if the end of the list is reached with no file having been appropriate (ie. found and overlapping) then an error is issued.

If, conversely, you wanted to ensure that at least the standard files were present, but that if all files were scanning without a hit a null transformation is applied you could use:

+nadgrids=conus,alaska,@ntv2_0.gsb,ntv1_can.dat,null

The null grid shift file is a valid grid shift file covering the whole world and applying no shift.

I still feel the default behavior (report an error if no covering grid shift file is found) is reasonable, but the behavior is fairly easily modified by changing the grid shift string. In the case of postgis, this could be done by replacing +datum=NAD27 with the above +nadgrids directive in the spatial_ref_sys table.

I'm not too keen on the patch approach.

[william....@…]

Thanks, Frank. I hadn't dug that deeply into proj (and thanks for all you do there and on GDAL). I can hardly do otherwise than gratefully embrace your analysis and conclusion. However …

1. It would be really nice if proj.4 error handling clearly distinguished between

missing grid-shift files and transformed points outside the range for which the grid shift is defined.

2. The 'second chance' logic in PostGIS's lwgeom_transform.c should be entirely

removed. It does not do what its comments claim, and, if it did, would conflict with the approach you recommend above.

3. It would be helpful to improve PostGIS warning/error messages as I suggest in the

next to last paragraph of my last message above, at least to the extent still relevant with your recommended approach.

I'll be happy to take a swipe at patches for items 2 and 3 above, though given my limited C coding skills, it would probably be better for the PostGIS development team to do so if resources are available.

— Bill

[mcayland]

I think that 2) and 3) should be done in time for PostGIS 1.4, so I'm marking this accordingly.

ATB,

Mark.

[mcayland]

William,

I've just committed the required changes for this to SVN trunk, including adding a new section to the manual to reflect this:

​http://postgis.refractions.net/documentation/manual-svn/ST_Transform.html

Since this is quite a major change in the way re-projection works, it's not really possible to backpatch to 1.3. However, if you wish to test then please download latest SVN trunk or if you wait a few days we hope to have a usable 1.4 beta available.

ATB,

Mark.

[mcayland]

<i>(No comment was entered for this change.)</i>

[william....@…]

​http://postgis.googlecode.com/issues/attachment?aid=141268354510428007&name=lwgeom_transform-1.3.3.diff

[william....@…]

Thanks a lot, Mark. I hadn't looked at recent updates of the PostGIS documentation.

The latest version is beautiful, from both an asethetic and functional perspective.

The last section of the updated ST_Transform page ('Configuring transformation behaviour') is all the detail one could ask for, even though most of the problematic action is at the proj.4 layer, below PostGIS itself. I'll try to give the SVN version a whirl this week and post results here.

— Bill

[william....@…]

​http://postgis.googlecode.com/issues/attachment?aid=-2669627169819104667&name=lwgeom_transform.c-1.3.3.diff

[mcayland]

William, are you happy with how current SVN trunk handles this?

ATB,

Mark.

[mcayland]

Assuming resolved due to lack of further feedback.