GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
geos.c
Go to the documentation of this file.
1 
16 #include <grass/config.h>
17 #include <grass/gis.h>
18 #include <grass/Vect.h>
19 #include <grass/glocale.h>
20 
21 #ifdef HAVE_GEOS
22 
23 static GEOSGeometry *Vect__read_line_geos(struct Map_info *, long, int *);
24 static GEOSCoordSequence *V1_read_line_geos(struct Map_info *, long, int *);
25 static GEOSCoordSequence *V2_read_line_geos(struct Map_info *, int);
26 static GEOSCoordSequence *read_polygon_points(struct Map_info *, int, int*);
27 
46 GEOSGeometry *Vect_read_line_geos(struct Map_info *Map, int line, int *type)
47 {
48  P_LINE *Line;
49 
50  G_debug(3, "Vect_read_line_geos(): line = %d", line);
51 
52  if (!VECT_OPEN(Map))
53  G_fatal_error("Vect_read_line_geos(): %s", _("vector map is not opened"));
54 
55  if (line < 1 || line > Map->plus.n_lines)
56  G_fatal_error(_("Vect_read_line_geos(): feature id %d is not reasonable "
57  "(max features in vector map <%s>: %d)"),
58  line, Vect_get_full_name(Map), Map->plus.n_lines);
59 
60  if (Map->format != GV_FORMAT_NATIVE)
61  G_fatal_error("Vect_read_line_geos(): %s", _("only native format supported"));
62 
63  Line = Map->plus.Line[line];
64  if (Line == NULL)
65  G_fatal_error("Vect_read_line_geos(): %s %d",
66  _("Attempt to read dead line"), line);
67 
68  return Vect__read_line_geos(Map, Line->offset, type);
69 }
70 
82 GEOSGeometry *Vect_read_area_geos(struct Map_info * Map, int area)
83 {
84  int i, nholes, isle;
85  GEOSGeometry *boundary, **holes;
86 
87  G_debug(3, "Vect_read_area_geos(): area = %d", area);
88 
89  boundary = GEOSGeom_createLinearRing(Vect_get_area_points_geos(Map, area));
90  if (!boundary) {
91  G_fatal_error(_("Vect_read_area_geos(): unable to read area id %d"),
92  area);
93  }
94 
95  nholes = Vect_get_area_num_isles(Map, area);
96  holes = (GEOSGeometry **) G_malloc(nholes * sizeof(GEOSGeometry *));
97  for (i = 0; i < nholes; i++) {
98  isle = Vect_get_area_isle(Map, area, i);
99  if (isle < 1) {
100  nholes--;
101  continue;
102  }
103  holes[i] = GEOSGeom_createLinearRing(Vect_get_isle_points_geos(Map, isle));
104  if (!(holes[i]))
105  G_fatal_error(_("Vect_read_area_geos(): unable to read isle id %d of area id %d"),
106  isle, area);
107  }
108 
109  return GEOSGeom_createPolygon(boundary, holes, nholes);
110 }
111 
128 GEOSGeometry *Vect_line_to_geos(struct Map_info *Map,
129  const struct line_pnts *points, int type)
130 {
131  int i, with_z;
132  GEOSGeometry *geom;
133  GEOSCoordSequence *pseq;
134 
135  G_debug(3, "Vect_line_to_geos(): type = %d", type);
136 
137  with_z = Vect_is_3d(Map);
138 
139  /* read only points / lines / boundaries */
140  if (!(type & (GV_POINT | GV_LINES)))
141  return NULL;
142 
143  if (type == GV_POINT) {
144  if (points->n_points != 1)
145  /* point is not valid */
146  return NULL;
147  }
148  else {
149  if (points->n_points < 2)
150  /* line/boundary is not valid */
151  return NULL;
152  }
153 
154  pseq = GEOSCoordSeq_create(points->n_points, with_z ? 3 : 2);
155 
156  for (i = 0; i < points->n_points; i++) {
157  GEOSCoordSeq_setX(pseq, i, points->x[i]);
158  GEOSCoordSeq_setY(pseq, i, points->y[i]);
159  if (with_z)
160  GEOSCoordSeq_setZ(pseq, i, points->z[i]);
161  }
162 
163  if (type == GV_POINT)
164  geom = GEOSGeom_createPoint(pseq);
165  else if (type == GV_LINE)
166  geom = GEOSGeom_createLineString(pseq);
167  else { /* boundary */
168  geom = GEOSGeom_createLineString(pseq);
169  if (GEOSisRing(geom)) {
170  /* GEOSGeom_destroy(geom); */
171  geom = GEOSGeom_createLinearRing(pseq);
172  }
173  }
174 
175  /* GEOSCoordSeq_destroy(pseq); */
176 
177  return geom;
178 }
179 
194 GEOSGeometry *Vect__read_line_geos(struct Map_info *Map, long offset, int *type)
195 {
196  int ftype;
197 
198  GEOSGeometry *geom;
199  GEOSCoordSequence *pseq;
200 
201  pseq = V1_read_line_geos(Map, offset, &ftype);
202  if (!pseq)
203  G_fatal_error(_("Unable to read line offset %ld"), offset);
204 
205  if (ftype & GV_POINT) {
206  G_debug(3, " geos_type = point");
207  geom = GEOSGeom_createPoint(pseq);
208  }
209  else if (ftype & GV_LINE) {
210  G_debug(3, " geos_type = linestring");
211  geom = GEOSGeom_createLineString(pseq);
212  }
213  else { /* boundary */
214  geom = GEOSGeom_createLineString(pseq);
215  if (GEOSisRing(geom)) {
216  /* GEOSGeom_destroy(geom); */
217  geom = GEOSGeom_createLinearRing(pseq);
218  G_debug(3, " geos_type = linearring");
219  }
220  else {
221  G_debug(3, " geos_type = linestring");
222  }
223  }
224 
225  /* GEOSCoordSeq_destroy(pseq); */
226 
227  if (type)
228  *type = ftype;
229 
230  return geom;
231 }
232 
245 GEOSCoordSequence *V2_read_line_geos(struct Map_info *Map, int line)
246 {
247  int ftype;
248  P_LINE *Line;
249 
250  G_debug(3, "V2_read_line_geos(): line = %d", line);
251 
252  Line = Map->plus.Line[line];
253 
254  if (Line == NULL)
255  G_fatal_error("V2_read_line_geos(): %s %d",
256  _("Attempt to read dead line"), line);
257 
258  return V1_read_line_geos(Map, Line->offset, &ftype);
259 }
260 
261 
278 GEOSCoordSequence *V1_read_line_geos(struct Map_info *Map, long offset, int *type)
279 {
280  int i, n_points;
281  int do_cats, n_cats;
282  char rhead, nc;
283  long size;
284  double *x, *y, *z;
285 
286  GEOSCoordSequence *pseq;
287 
288  G_debug(3, "V1_read_line_geos(): offset = %ld", offset);
289 
290  Map->head.last_offset = offset;
291 
292  /* reads must set in_head, but writes use default */
293  dig_set_cur_port(&(Map->head.port));
294 
295  dig_fseek(&(Map->dig_fp), offset, 0);
296 
297  if (0 >= dig__fread_port_C(&rhead, 1, &(Map->dig_fp)))
298  return NULL; /* end of file */
299 
300  if (!(rhead & 0x01)) /* dead line */
301  return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
302 
303  if (rhead & 0x02) /* categories exists */
304  do_cats = 1; /* do not return here let file offset moves forward to next */
305  else /* line */
306  do_cats = 0;
307 
308  rhead >>= 2;
309  *type = dig_type_from_store((int) rhead);
310 
311  /* read only points / lines / boundaries */
312  if (!(*type & (GV_POINT | GV_LINES)))
313  return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
314 
315  /* skip categories */
316  if (do_cats) {
317  if (Map->head.Version_Minor == 1) { /* coor format 5.1 */
318  if (0 >= dig__fread_port_I(&n_cats, 1, &(Map->dig_fp)))
319  return NULL;
320  }
321  else { /* coor format 5.0 */
322  if (0 >= dig__fread_port_C(&nc, 1, &(Map->dig_fp)))
323  return NULL;
324  n_cats = (int) nc;
325  }
326  G_debug(3, " n_cats = %d", n_cats);
327 
328  if (Map->head.Version_Minor == 1) { /* coor format 5.1 */
329  size = (2 * PORT_INT) * n_cats;
330  }
331  else { /* coor format 5.0 */
332  size = (PORT_SHORT + PORT_INT) * n_cats;
333  }
334  dig_fseek(&(Map->dig_fp), size, SEEK_CUR);
335  }
336 
337  if (*type & GV_POINTS) {
338  n_points = 1;
339  }
340  else {
341  if (0 >= dig__fread_port_I(&n_points, 1, &(Map->dig_fp)))
342  return NULL;
343  }
344 
345  G_debug(3, " n_points = %d dim = %d", n_points, (Map->head.with_z) ? 3 : 2);
346 
347  pseq = GEOSCoordSeq_create(n_points, (Map->head.with_z) ? 3 : 2);
348 
349  x = (double *) G_malloc(n_points * sizeof(double));
350  y = (double *) G_malloc(n_points * sizeof(double));
351  if (Map->head.with_z)
352  z = (double *) G_malloc(n_points * sizeof(double));
353  else
354  z = NULL;
355 
356  if (0 >= dig__fread_port_D(x, n_points, &(Map->dig_fp)))
357  return NULL; /* end of file */
358 
359  if (0 >= dig__fread_port_D(y, n_points, &(Map->dig_fp)))
360  return NULL; /* end of file */
361 
362  if (Map->head.with_z) {
363  if (0 >= dig__fread_port_D(z, n_points, &(Map->dig_fp)))
364  return NULL; /* end of file */
365 
366  }
367 
368  for (i = 0; i < n_points; i++) {
369  GEOSCoordSeq_setX(pseq, i, x[i]);
370  GEOSCoordSeq_setY(pseq, i, y[i]);
371  if (Map->head.with_z)
372  GEOSCoordSeq_setZ(pseq, i, z[i]);
373  }
374 
375  G_debug(3, " off = %ld", dig_ftell(&(Map->dig_fp)));
376 
377  G_free((void *) x);
378  G_free((void *) y);
379  if (z)
380  G_free((void *) z);
381 
382  return pseq;
383 }
384 
399 GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *Map, int area)
400 {
401  struct Plus_head *Plus;
402  P_AREA *Area;
403 
404  G_debug(3, "Vect_get_area_points_geos(): area = %d", area);
405 
406  Plus = &(Map->plus);
407  Area = Plus->Area[area];
408 
409  if (Area == NULL) { /* dead area */
410  G_warning(_("Attempt to read points of nonexistent area id %d"), area);
411  return NULL; /* error , because we should not read dead areas */
412  }
413 
414  return read_polygon_points(Map, Area->n_lines, Area->lines);
415 }
416 
430 GEOSCoordSequence *Vect_get_isle_points_geos(struct Map_info *Map, int isle)
431 {
432  struct Plus_head *Plus;
433  P_ISLE *Isle;
434 
435  G_debug(3, "Vect_get_isle_points_geos(): isle = %d", isle);
436 
437  Plus = &(Map->plus);
438  Isle = Plus->Isle[isle];
439 
440  return read_polygon_points(Map, Isle->n_lines, Isle->lines);
441 }
442 
443 GEOSCoordSequence *read_polygon_points(struct Map_info *Map, int n_lines, int *lines)
444 {
445  int i, j, k;
446  int line, aline;
447  unsigned int n_points, n_points_shell;
448  double x, y, z;
449  int *dir;
450 
451  GEOSCoordSequence **pseq, *pseq_shell;
452 
453  G_debug(3, " n_lines = %d", n_lines);
454  pseq = (GEOSCoordSequence **) G_malloc(n_lines * sizeof(GEOSCoordSequence *));
455  dir = (int*) G_malloc(n_lines * sizeof(int));
456 
457  n_points_shell = 0;
458  for (i = 0; i < n_lines; i++) {
459  line = lines[i];
460  aline = abs(line);
461  G_debug(3, " append line(%d) = %d", i, line);
462 
463  if (line > 0)
464  dir[i] = GV_FORWARD;
465  else
466  dir[i] = GV_BACKWARD;
467 
468  pseq[i] = V2_read_line_geos(Map, aline);
469  if (!(pseq[i])) {
470  G_fatal_error(_("Unable to read feature id %d"), aline);
471  }
472 
473  GEOSCoordSeq_getSize(pseq[i], &n_points);
474  G_debug(3, " line n_points = %d", n_points);
475  n_points_shell += n_points;
476  }
477 
478  /* create shell (outer ring) */
479  pseq_shell = GEOSCoordSeq_create(n_points_shell, Map->head.with_z ? 3 : 2);
480  k = 0;
481  for (i = 0; i < n_lines; i++) {
482  GEOSCoordSeq_getSize(pseq[i], &n_points);
483  if (dir[i] == GV_FORWARD) {
484  for (j = 0; j < (int) n_points; j++, k++) {
485  GEOSCoordSeq_getX(pseq[i], j, &x);
486  GEOSCoordSeq_setX(pseq_shell, k, x);
487 
488  GEOSCoordSeq_getY(pseq[i], j, &y);
489  GEOSCoordSeq_setY(pseq_shell, k, y);
490 
491  if (Map->head.with_z) {
492  GEOSCoordSeq_getY(pseq[i], j, &z);
493  GEOSCoordSeq_setZ(pseq_shell, k, z);
494  }
495  }
496  }
497  else { /* GV_BACKWARD */
498  for (j = (int) n_points - 1; j > -1; j--, k++) {
499  GEOSCoordSeq_getX(pseq[i], j, &x);
500  GEOSCoordSeq_setX(pseq_shell, k, x);
501 
502  GEOSCoordSeq_getY(pseq[i], j, &y);
503  GEOSCoordSeq_setY(pseq_shell, k, y);
504 
505  if (Map->head.with_z) {
506  GEOSCoordSeq_getY(pseq[i], j, &z);
507  GEOSCoordSeq_setZ(pseq_shell, k, z);
508  }
509  }
510  }
511  }
512 
513  G_free((void *) pseq);
514  G_free((void *) dir);
515 
516  return pseq_shell;
517 }
518 #endif /* HAVE_GEOS */
GEOSCoordSequence * Vect_get_area_points_geos(struct Map_info *Map, int area)
Returns the polygon array of points, i.e. outer ring (shell)
Definition: geos.c:399
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
int dig_set_cur_port(struct Port_info *port)
Definition: portable.c:640
GEOSCoordSequence * Vect_get_isle_points_geos(struct Map_info *Map, int isle)
Returns the polygon (isle) array of points (inner ring)
Definition: geos.c:430
int Vect_get_area_num_isles(struct Map_info *Map, int area)
Returns number of isles for area.
int Vect_get_area_isle(struct Map_info *Map, int area, int isle)
Returns isle for area.
int y
Definition: plot.c:34
GEOSGeometry * Vect_read_line_geos(struct Map_info *Map, int line, int *type)
Read vector feature and stores it as GEOSGeometry instance.
Definition: geos.c:46
const char * Vect_get_full_name(struct Map_info *Map)
Get full map name.
tuple size
value.Bind(wx.EVT_TEXT, self.OnVolumeIsosurfMap)
Definition: tools.py:2334
int dig_fseek(GVFILE *file, long offset, int whence)
Set GVFILE position.
Definition: file.c:60
int Vect_is_3d(struct Map_info *Map)
Check if vector map is 3D (with z)
int GV_LINES
Definition: vdigit/main.py:24
int dig_type_from_store(int stype)
Convert type from store type.
int
Definition: g3dcolor.c:48
int dig__fread_port_D(double *buf, int cnt, GVFILE *fp)
Definition: portable.c:75
GEOSGeometry * Vect_line_to_geos(struct Map_info *Map, const struct line_pnts *points, int type)
Create GEOSGeometry of given type from feature points.
Definition: geos.c:128
int dig__fread_port_C(char *buf, int cnt, GVFILE *fp)
Definition: portable.c:347
return NULL
Definition: dbfopen.c:1394
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
tuple Map
Definition: render.py:1310
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
GEOSGeometry * Vect_read_area_geos(struct Map_info *Map, int area)
Read vector area and stores it as GEOSGeometry instance (polygon)
Definition: geos.c:82
int G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
long dig_ftell(GVFILE *file)
Get GVFILE position.
Definition: file.c:36
int dig__fread_port_I(int *buf, int cnt, GVFILE *fp)
Definition: portable.c:207