PostGIS  2.5.0dev-r@@SVN_REVISION@@
liblwgeom/lwgeom_transform.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright (C) 2001-2003 Refractions Research Inc.
22  *
23  **********************************************************************/
24 
25 
26 #include "../postgis_config.h"
27 #include "liblwgeom_internal.h"
28 #include "lwgeom_log.h"
29 #include <string.h>
30 
31 
33 static void
35 {
36  pt->x *= M_PI/180.0;
37  pt->y *= M_PI/180.0;
38 }
39 
41 static void
43 {
44  pt->x *= 180.0/M_PI;
45  pt->y *= 180.0/M_PI;
46 }
47 
52 int
53 ptarray_transform(POINTARRAY *pa, projPJ inpj, projPJ outpj)
54 {
55  uint32_t i;
56  POINT4D p;
57 
58  for ( i = 0; i < pa->npoints; i++ )
59  {
60  getPoint4d_p(pa, i, &p);
61  if ( ! point4d_transform(&p, inpj, outpj) ) return LW_FAILURE;
62  ptarray_set_point4d(pa, i, &p);
63  }
64 
65  return LW_SUCCESS;
66 }
67 
68 
73 int
74 lwgeom_transform(LWGEOM *geom, projPJ inpj, projPJ outpj)
75 {
76  uint32_t i;
77 
78  /* No points to transform in an empty! */
79  if ( lwgeom_is_empty(geom) )
80  return LW_SUCCESS;
81 
82  switch(geom->type)
83  {
84  case POINTTYPE:
85  case LINETYPE:
86  case CIRCSTRINGTYPE:
87  case TRIANGLETYPE:
88  {
89  LWLINE *g = (LWLINE*)geom;
90  if ( ! ptarray_transform(g->points, inpj, outpj) ) return LW_FAILURE;
91  break;
92  }
93  case POLYGONTYPE:
94  {
95  LWPOLY *g = (LWPOLY*)geom;
96  for ( i = 0; i < g->nrings; i++ )
97  {
98  if ( ! ptarray_transform(g->rings[i], inpj, outpj) ) return LW_FAILURE;
99  }
100  break;
101  }
102  case MULTIPOINTTYPE:
103  case MULTILINETYPE:
104  case MULTIPOLYGONTYPE:
105  case COLLECTIONTYPE:
106  case COMPOUNDTYPE:
107  case CURVEPOLYTYPE:
108  case MULTICURVETYPE:
109  case MULTISURFACETYPE:
111  case TINTYPE:
112  {
113  LWCOLLECTION *g = (LWCOLLECTION*)geom;
114  for ( i = 0; i < g->ngeoms; i++ )
115  {
116  if ( ! lwgeom_transform(g->geoms[i], inpj, outpj) ) return LW_FAILURE;
117  }
118  break;
119  }
120  default:
121  {
122  lwerror("lwgeom_transform: Cannot handle type '%s'",
123  lwtype_name(geom->type));
124  return LW_FAILURE;
125  }
126  }
127  return LW_SUCCESS;
128 }
129 
130 int
131 point4d_transform(POINT4D *pt, projPJ srcpj, projPJ dstpj)
132 {
133  int* pj_errno_ref;
134  POINT4D orig_pt;
135 
136  /* Make a copy of the input point so we can report the original should an error occur */
137  orig_pt.x = pt->x;
138  orig_pt.y = pt->y;
139  orig_pt.z = pt->z;
140 
141  if (pj_is_latlong(srcpj)) to_rad(pt) ;
142 
143  LWDEBUGF(4, "transforming POINT(%f %f) from '%s' to '%s'", orig_pt.x, orig_pt.y, pj_get_def(srcpj,0), pj_get_def(dstpj,0));
144 
145  /* Perform the transform */
146  pj_transform(srcpj, dstpj, 1, 0, &(pt->x), &(pt->y), &(pt->z));
147 
148  /* For NAD grid-shift errors, display an error message with an additional hint */
149  pj_errno_ref = pj_get_errno_ref();
150 
151  if (*pj_errno_ref != 0)
152  {
153  if (*pj_errno_ref == -38)
154  {
155  lwnotice("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.");
156  lwerror("transform: couldn't project point (%g %g %g): %s (%d)",
157  orig_pt.x, orig_pt.y, orig_pt.z, pj_strerrno(*pj_errno_ref), *pj_errno_ref);
158  return 0;
159  }
160  else
161  {
162  lwerror("transform: couldn't project point (%g %g %g): %s (%d)",
163  orig_pt.x, orig_pt.y, orig_pt.z, pj_strerrno(*pj_errno_ref), *pj_errno_ref);
164  return 0;
165  }
166  }
167 
168  if (pj_is_latlong(dstpj)) to_dec(pt);
169  return 1;
170 }
171 
172 projPJ
173 lwproj_from_string(const char *str1)
174 {
175  int t;
176  char *params[1024]; /* one for each parameter */
177  char *loc;
178  char *str;
179  size_t slen;
180  projPJ result;
181 
182 
183  if (str1 == NULL) return NULL;
184 
185  slen = strlen(str1);
186 
187  if (slen == 0) return NULL;
188 
189  str = lwalloc(slen+1);
190  strcpy(str, str1);
191 
192  /*
193  * first we split the string into a bunch of smaller strings,
194  * based on the " " separator
195  */
196 
197  params[0] = str; /* 1st param, we'll null terminate at the " " soon */
198 
199  loc = str;
200  t = 1;
201  while ((loc != NULL) && (*loc != 0) )
202  {
203  loc = strchr(loc, ' ');
204  if (loc != NULL)
205  {
206  *loc = 0; /* null terminate */
207  params[t] = loc+1;
208  loc++; /* next char */
209  t++; /*next param */
210  }
211  }
212 
213  if (!(result=pj_init(t, params)))
214  {
215  lwfree(str);
216  return NULL;
217  }
218  lwfree(str);
219  return result;
220 }
221 
222 
223 
double x
Definition: liblwgeom.h:351
#define LINETYPE
Definition: liblwgeom.h:85
static void to_dec(POINT4D *pt)
convert radians to decimal degress
#define MULTICURVETYPE
Definition: liblwgeom.h:94
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
void lwfree(void *mem)
Definition: lwutil.c:244
int lwgeom_transform(LWGEOM *geom, projPJ inpj, projPJ outpj)
Transform given SERIALIZED geometry from inpj projection to outpj projection.
#define POLYGONTYPE
Definition: liblwgeom.h:86
#define CURVEPOLYTYPE
Definition: liblwgeom.h:93
#define COMPOUNDTYPE
Definition: liblwgeom.h:92
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
#define LW_SUCCESS
Definition: liblwgeom.h:79
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:425
#define TRIANGLETYPE
Definition: liblwgeom.h:97
#define POLYHEDRALSURFACETYPE
Definition: liblwgeom.h:96
uint32_t ngeoms
Definition: liblwgeom.h:506
uint32_t nrings
Definition: liblwgeom.h:454
#define LW_FAILURE
Definition: liblwgeom.h:78
unsigned int uint32_t
Definition: uthash.h:78
projPJ lwproj_from_string(const char *str1)
Get a projection from a string representation.
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
int point4d_transform(POINT4D *pt, projPJ srcpj, projPJ dstpj)
LWGEOM ** geoms
Definition: liblwgeom.h:508
#define TINTYPE
Definition: liblwgeom.h:98
POINTARRAY ** rings
Definition: liblwgeom.h:456
double z
Definition: liblwgeom.h:351
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:89
#define MULTISURFACETYPE
Definition: liblwgeom.h:95
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
uint8_t type
Definition: liblwgeom.h:395
#define CIRCSTRINGTYPE
Definition: liblwgeom.h:91
void * lwalloc(size_t size)
Definition: lwutil.c:229
static void to_rad(POINT4D *pt)
convert decimal degress to radians
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1386
double y
Definition: liblwgeom.h:351
#define MULTILINETYPE
Definition: liblwgeom.h:88
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
int ptarray_transform(POINTARRAY *pa, projPJ inpj, projPJ outpj)
Transform given POINTARRAY from inpj projection to outpj projection.
POINTARRAY * points
Definition: liblwgeom.h:421
uint32_t npoints
Definition: liblwgeom.h:370