self-hosted/postgis
self-hosted/postgis
1
0
Fork 0
mirror of https://git.osgeo.org/gitea/postgis/postgis synced 2024-10-25 17:42:38 +00:00
Code Issues Packages Projects Releases Wiki Activity

Fix for GBT#89: transform() grid-shift 2nd chance logic defective. Remove the 2nd chance logic completely and allow the user to configure the behaviour using the standard PROJ.4 +nadgrids parameter. I've added a section to the ST_Transform() section of the manual which gives an example of how you can do this.

Browse source 
git-svn-id: http://svn.osgeo.org/postgis/trunk@3825 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Mark Cave-Ayland 2009-03-10 15:06:42 +00:00
parent 9b858e03dc
commit 8ef685afd9
2 changed files with 39 additions and 81 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

12
doc/reference_new.xml
View file

@ -6765,6 +6765,18 @@ CREATE INDEX idx_the_geom_26986_parcels
</programlisting>
</refsection>
<refsection>
<title>Configuring transformation behaviour</title>
<para>Sometimes coordinate transformation involving a grid-shift can fail, for example if PROJ.4 has not been built with grid-shift files or the coordinate does not lie within the range for which the grid shift is defined. By default, PostGIS will throw an error if a grid shift file is not present, but this behaviour can be configured on a per-SRID basis by altering the proj4text value within the spatial_ref_sys table.</para>
<para>For example, the proj4text parameter +datum=NAD87 is a shorthand form for the following +nadgrids parameter:</para>
<programlisting>+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat</programlisting>
<para>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.</para>
<para>If, conversely, you wanted to ensure that at least the standard files were present, but that if all files were scanned without a hit a null transformation is applied you could use:</para>
<programlisting>+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat,null</programlisting>
<para>The null grid shift file is a valid grid shift file covering the whole world and applying no shift. So for a complete example, if you wanted to alter PostGIS so that transformations to SRID 4267 that didn't lie within the correct range did not throw an ERROR, you would use the following:</para>
<programlisting>UPDATE spatial_ref_sys SET proj4text = '+proj=longlat +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat,null +no_defs' WHERE srid = 4267;</programlisting>
</refsection>
<!-- Optionally add a "See Also" section -->
<refsection>
<title>See Also</title>

108
postgis/lwgeom_transform.c
View file

@ -521,69 +521,6 @@ void SetPROJ4LibPath(void)
}
/*
* This is *exactly* the same as PROJ.4's pj_transform(),
* but it doesn't do the datum shift.
*/
int
pj_transform_nodatum(PJ *srcdefn, PJ *dstdefn, long point_count,
int point_offset, double *x, double *y, double *z )
{
long i;
/*int need_datum_shift; */
int* pj_errno_ref;
if( point_offset == 0 )
point_offset = 1;
if( !srcdefn->is_latlong )
{
for( i = 0; i < point_count; i++ )
{
XY projected_loc;
LP geodetic_loc;
projected_loc.u = x[point_offset*i];
projected_loc.v = y[point_offset*i];
geodetic_loc = pj_inv( projected_loc, srcdefn );
pj_errno_ref = pj_get_errno_ref();
if( (*pj_errno_ref) != 0 )
return *pj_errno_ref;
x[point_offset*i] = geodetic_loc.u;
y[point_offset*i] = geodetic_loc.v;
}
}
if( !dstdefn->is_latlong )
{
for( i = 0; i < point_count; i++ )
{
XY projected_loc;
LP geodetic_loc;
geodetic_loc.u = x[point_offset*i];
geodetic_loc.v = y[point_offset*i];
projected_loc = pj_fwd( geodetic_loc, dstdefn );
pj_errno_ref = pj_get_errno_ref();
if( (*pj_errno_ref) != 0 )
return *pj_errno_ref;
x[point_offset*i] = projected_loc.u;
y[point_offset*i] = projected_loc.v;
}
}
return 0;
}
/* convert decimal degress to radians */
void
to_rad(POINT4D *pt)
@ -1006,27 +943,36 @@ int
transform_point(POINT4D *pt, PJ *srcpj, PJ *dstpj)
{
int* pj_errno_ref;
POINT4D orig_pt;
if (srcpj->is_latlong) to_rad(pt);
pj_transform(srcpj, dstpj, 1, 2, &(pt->x), &(pt->y), &(pt->z));
pj_errno_ref = pj_get_errno_ref();
if (*pj_errno_ref)
{
if ((*pj_errno_ref) == -38) /*2nd chance */
{
elog(WARNING, "transform: %i (%s)",
*pj_errno_ref, pj_strerrno(*pj_errno_ref));
/*couldnt do nadshift - do it without the datum */
pj_transform_nodatum(srcpj, dstpj, 1, 2,
&(pt->x), &(pt->y), NULL);
}
/* Make a copy of the input point so we can report the original should an error occur */
orig_pt.x = pt->x;
orig_pt.y = pt->y;
orig_pt.z = pt->z;
pj_errno_ref = pj_get_errno_ref();
if ((*pj_errno_ref))
if (srcpj->is_latlong) to_rad(pt);
/* Perform the transform */
pj_transform(srcpj, dstpj, 1, 0, &(pt->x), &(pt->y), &(pt->z));
/* For NAD grid-shift errors, display an error message with an additional hint */
pj_errno_ref = pj_get_errno_ref();
if (*pj_errno_ref != 0)
{
if (*pj_errno_ref == -38)
{
elog(ERROR,"transform: couldn't project point: %i (%s)",
*pj_errno_ref, pj_strerrno(*pj_errno_ref));
ereport(ERROR, (
errmsg_internal("transform: couldn't project point (%g %g %g): %s (%d)",
orig_pt.x, orig_pt.y, orig_pt.z, pj_strerrno(*pj_errno_ref), *pj_errno_ref),
errhint("PostGIS was unable to transform the point because either no grid shift files were found, or the point does not lie within the range for which the grid shift is defined. Refer to the ST_Transform() section of the PostGIS manual for details on how to configure PostGIS to alter this behaviour.")
));
return 0;
}
else
{
elog(ERROR, "transform: couldn't project point (%g %g %g): %s (%d)",
orig_pt.x, orig_pt.y, orig_pt.z, pj_strerrno(*pj_errno_ref), *pj_errno_ref);
return 0;
}
}