boundary splitting distance: account for multipolygons

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@43040 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/v.in.ogr/geom.c

@@ -158,6 +158,37 @@ centroid(OGRGeometryH hGeom, CENTR * Centr, struct spatial_index * Sindex, int f
     return 0;
 }
 
+/* count polygons and isles */
+int
+poly_count(OGRGeometryH hGeom)
+{
+    int i, nr, ret;
+    OGRwkbGeometryType eType;
+    OGRGeometryH hRing;
+
+    eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
+
+    if (eType == wkbPolygon) {
+	G_debug(3, "Polygon");
+	nr = OGR_G_GetGeometryCount(hGeom);
+	n_polygon_boundaries += nr;
+
+    }
+    else if (eType == wkbGeometryCollection || eType == wkbMultiPolygon) {
+	G_debug(3, "GeometryCollection or MultiPolygon");
+	nr = OGR_G_GetGeometryCount(hGeom);
+	for (i = 0; i < nr; i++) {
+	    hRing = OGR_G_GetGeometryRef(hGeom, i);
+
+	    ret = poly_count(hRing);
+	    if (ret == -1) {
+		G_warning(_("Cannot read part of geometry"));
+	    }
+	}
+    }
+    return 0;
+}
+
 /* Write geometry to output map */
 int
 geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,

vector/v.in.ogr/global.h

@@ -25,6 +25,7 @@
 
 
 extern int n_polygons;
+extern int n_polygon_boundaries;
 extern double split_distance;
 
 /* centroid structure */

vector/v.in.ogr/main.c

@@ -38,12 +38,14 @@
 #endif
 
 int n_polygons;
+int n_polygon_boundaries;
 double split_distance;
 
 int geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,
 	 double min_area, int type, int mk_centr);
 int centroid(OGRGeometryH hGeom, CENTR * Centr, struct spatial_index * Sindex,
 	     int field, int cat, double min_area, int type);
+int poly_count(OGRGeometryH hGeom);
 
 int main(int argc, char *argv[])
 {
@@ -820,16 +822,38 @@ int main(int argc, char *argv[])
 	OGR_L_ResetReading(Ogr_layer);
 	unsigned int n_features = 0, feature_count = 0;
 
+	n_polygon_boundaries = 0;
 	n_features = OGR_L_GetFeatureCount(Ogr_layer, 1);
-	if (split_distance > -0.5 && n_features > 500) {
+
+	/* estimate distance for boundary splitting --> */
+
+	/* count polygons and isles */
+	G_message(_("Counting polygons for %d features..."), n_features);
+	while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
+	    G_percent(feature_count++, n_features, 1);	/* show something happens */
+	    /* Geometry */
+	    Ogr_geometry = OGR_F_GetGeometryRef(Ogr_feature);
+	    if (Ogr_geometry != NULL) {
+		poly_count(Ogr_geometry);
+	    }
+	    OGR_F_Destroy(Ogr_feature);
+	}
+	/* rewind layer */
+	OGR_L_ResetReading(Ogr_layer);
+	feature_count = 0;
+
+	G_debug(0, "n polygon boundaries: %d", n_polygon_boundaries);
+	if (split_distance > -0.5 && n_polygon_boundaries > 50) {
 	    split_distance =
-		area_size / log(OGR_L_GetFeatureCount(Ogr_layer, 1));
+		area_size / log(n_features);
 	    /* divisor is the handle: increase divisor to decrease split_distance */
 	    split_distance = split_distance / 4.;
 	    G_debug(1, "root of area size: %f", area_size);
 	    G_verbose_message(_("Boundary splitting distance in map units: %G"),
 		      split_distance);
 	}
+	/* <-- estimate distance for boundary splitting */
+	
 	G_important_message(_("Importing map %d features..."), n_features);
 	while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
 	    G_percent(feature_count++, n_features, 1);	/* show something happens */