← Back to branch summary

~ubuntu-branches/ubuntu/precise/grass/precise

« back to all changes in this revision

Viewing changes to vector/v.in.ogr/main.c

  • Committer: Bazaar Package Importer
  • Author(s): Francesco Paolo Lovergine
  • Date: 2011-04-13 17:08:41 UTC
  • mfrom: (8.1.7 sid)
  • Revision ID: james.westby@ubuntu.com-20110413170841-ss1t9bic0d0uq0gz
Tags: 6.4.1-1
* New upstream version.
* Now build-dep on libjpeg-dev and current libreadline6-dev.
* Removed patch swig: obsolete.
* Policy bumped to 3.9.2, without changes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
vector/v.in.ogr/main.c
41
41
         double min_area, int type, int mk_centr);
42
42
int centroid(OGRGeometryH hGeom, CENTR * Centr, SPATIAL_INDEX * Sindex,
43
43
             int field, int cat, double min_area, int type);
 
44
int poly_count(OGRGeometryH hGeom);
44
45
 
45
46
int main(int argc, char *argv[])
46
47
{
55
56
    struct Flag *list_flag, *no_clean_flag, *z_flag, *notab_flag,
56
57
        *region_flag;
57
58
    struct Flag *over_flag, *extend_flag, *formats_flag, *tolower_flag;
58
 
    char buf[2000], namebuf[2000];
 
59
    char buf[2000], namebuf[2000], tempvect[2000];
59
60
    char *separator;
60
61
    struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
61
62
    struct Key_Value *proj_info, *proj_units;
63
64
    char error_msg[8192];
64
65
 
65
66
    /* Vector */
66
 
    struct Map_info Map;
 
67
    struct Map_info Map, Tmp;
67
68
    int cat;
68
69
 
69
70
    /* Attributes */
92
93
    int navailable_layers;
93
94
    int layer_id;
94
95
    int overwrite;
 
96
    double area_size = 0.;
95
97
 
96
98
    G_gisinit(argv[0]);
97
99
 
456
458
        cellhd.tb_res = 1.;
457
459
    }
458
460
 
 
461
    /* split boundaries only when cleaning */
 
462
    if (no_clean_flag->answer) {
 
463
        split_distance = -1.;
 
464
    }
 
465
    else {
 
466
        split_distance = 0.;
 
467
        area_size =
 
468
            sqrt((cellhd.east - cellhd.west) * (cellhd.north - cellhd.south));
 
469
    }
 
470
 
459
471
    /* Fetch input map projection in GRASS form. */
460
472
    proj_info = NULL;
461
473
    proj_units = NULL;
586
598
    db_init_string(&strval);
587
599
 
588
600
    /* open output vector */
589
 
    if (z_flag->answer)
 
601
    /* open temporary vector, do the work in the temporary vector
 
602
     * at the end copy alive lines to output vector
 
603
     * in case of polygons this reduces the coor file size by a factor of 2 to 5
 
604
     * only needed for polygons, but the presence of polygons can be detected
 
605
     * only during OGR feature import, not before */
 
606
    sprintf(buf, "%s", out_opt->answer);
 
607
    /* strip any @mapset from vector output name */
 
608
    G_find_vector(buf, G_mapset());
 
609
    sprintf(tempvect, "%s_tmp", buf);
 
610
    G_verbose_message(_("Using temporary vector <%s>"), tempvect);
 
611
    if (z_flag->answer) {
590
612
        Vect_open_new(&Map, out_opt->answer, 1);
591
 
    else
 
613
        Vect_open_new(&Tmp, tempvect, 1);
 
614
    }
 
615
    else {
592
616
        Vect_open_new(&Map, out_opt->answer, 0);
 
617
        Vect_open_new(&Tmp, tempvect, 0);
 
618
    }
593
619
 
594
620
    Vect_hist_command(&Map);
595
621
 
622
648
            if (ncnames > 0) {
623
649
                cat_col_name = cnames_opt->answers[0];
624
650
            }
625
 
            Vect_map_add_dblink(&Map, layer + 1, NULL, Fi->table,
 
651
            Vect_map_add_dblink(&Map, layer + 1, layer_names[layer], Fi->table,
626
652
                                cat_col_name, Fi->database, Fi->driver);
627
653
 
628
654
            ncols = OGR_FD_GetFieldCount(Ogr_featuredefn);
782
808
        cat = 1;
783
809
        nogeom = 0;
784
810
        OGR_L_ResetReading(Ogr_layer);
785
 
        G_important_message(_("Importing map %d features..."),
786
 
                            OGR_L_GetFeatureCount(Ogr_layer, 1));
 
811
        unsigned int n_features = 0, feature_count = 0;
 
812
 
 
813
        n_polygon_boundaries = 0;
 
814
        n_features = OGR_L_GetFeatureCount(Ogr_layer, 1);
 
815
 
 
816
        /* estimate distance for boundary splitting --> */
 
817
 
 
818
        if (split_distance > -0.5) {
 
819
            /* count polygons and isles */
 
820
            G_message(_("Counting polygons for %d features..."), n_features);
 
821
            while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
 
822
                G_percent(feature_count++, n_features, 1);      /* show something happens */
 
823
                /* Geometry */
 
824
                Ogr_geometry = OGR_F_GetGeometryRef(Ogr_feature);
 
825
                if (Ogr_geometry != NULL) {
 
826
                    poly_count(Ogr_geometry);
 
827
                }
 
828
                OGR_F_Destroy(Ogr_feature);
 
829
            }
 
830
            /* rewind layer */
 
831
            OGR_L_ResetReading(Ogr_layer);
 
832
            feature_count = 0;
 
833
        }
 
834
 
 
835
        G_debug(1, "n polygon boundaries: %d", n_polygon_boundaries);
 
836
        if (split_distance > -0.5 && n_polygon_boundaries > 50) {
 
837
            split_distance =
 
838
                area_size / log(n_polygon_boundaries);
 
839
            /* divisor is the handle: increase divisor to decrease split_distance */
 
840
            split_distance = split_distance / 4.;
 
841
            G_debug(1, "root of area size: %f", area_size);
 
842
            G_verbose_message(_("Boundary splitting distance in map units: %G"),
 
843
                      split_distance);
 
844
        }
 
845
        /* <-- estimate distance for boundary splitting */
 
846
        
 
847
        G_important_message(_("Importing map %d features..."), n_features);
787
848
        while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
 
849
            G_percent(feature_count++, n_features, 1);  /* show something happens */
788
850
            /* Geometry */
789
851
            Ogr_geometry = OGR_F_GetGeometryRef(Ogr_feature);
790
852
            if (Ogr_geometry == NULL) {
795
857
                if (dim > 2)
796
858
                    with_z = 1;
797
859
 
798
 
                geom(Ogr_geometry, &Map, layer + 1, cat, min_area, type,
 
860
                geom(Ogr_geometry, &Tmp, layer + 1, cat, min_area, type,
799
861
                     no_clean_flag->answer);
800
862
            }
801
863
 
869
931
            OGR_F_Destroy(Ogr_feature);
870
932
            cat++;
871
933
        }
 
934
        G_percent(n_features, n_features, 1);   /* finish it */
872
935
 
873
936
        if (!notab_flag->answer) {
874
937
            db_commit_transaction(driver);
885
948
    G_message("%s", separator);
886
949
 
887
950
    /* TODO: is it necessary to build here? probably not, consumes time */
888
 
    Vect_build(&Map);
 
951
    /* GV_BUILD_BASE is sufficient to toggle boundary cleaning */
 
952
    Vect_build_partial(&Tmp, GV_BUILD_BASE);
889
953
 
890
954
    if (!no_clean_flag->answer &&
891
 
        Vect_get_num_primitives(&Map, GV_BOUNDARY) > 0) {
 
955
        Vect_get_num_primitives(&Tmp, GV_BOUNDARY) > 0) {
892
956
        int ret, centr, ncentr, otype, n_overlaps, n_nocat;
893
957
        CENTR *Centr;
894
958
        SPATIAL_INDEX si;
902
966
        G_message("%s", separator);
903
967
        G_warning(_("Cleaning polygons, result is not guaranteed!"));
904
968
 
905
 
        Vect_set_release_support(&Map);
906
 
        Vect_close(&Map);
907
 
        Vect_open_update(&Map, out_opt->answer, G_mapset());
908
 
        Vect_build_partial(&Map, GV_BUILD_BASE);        /* Downgrade topo */
909
 
 
910
969
        if (snap >= 0) {
911
970
            G_message("%s", separator);
912
971
            G_message(_("Snap boundaries (threshold = %.3e):"), snap);
913
 
            Vect_snap_lines(&Map, GV_BOUNDARY, snap, NULL);
 
972
            Vect_snap_lines(&Tmp, GV_BOUNDARY, snap, NULL);
914
973
        }
915
974
 
916
975
        /* It is not to clean to snap centroids, but I have seen data with 2 duplicate polygons
924
983
 
925
984
        G_message("%s", separator);
926
985
        G_message(_("Break polygons:"));
927
 
        Vect_break_polygons(&Map, GV_BOUNDARY, NULL);
 
986
        Vect_break_polygons(&Tmp, GV_BOUNDARY, NULL);
928
987
 
929
 
        /* It is important to remove also duplicate centroids in case of duplicate imput polygons */
 
988
        /* It is important to remove also duplicate centroids in case of duplicate input polygons */
930
989
        G_message("%s", separator);
931
990
        G_message(_("Remove duplicates:"));
932
 
        Vect_remove_duplicates(&Map, GV_BOUNDARY | GV_CENTROID, NULL);
 
991
        Vect_remove_duplicates(&Tmp, GV_BOUNDARY | GV_CENTROID, NULL);
 
992
 
 
993
        /* split boundaries here ? */
933
994
 
934
995
        /* Vect_clean_small_angles_at_nodes() can change the geometry so that new intersections
935
996
         * are created. We must call Vect_break_lines(), Vect_remove_duplicates()
936
 
         * and Vect_clean_small_angles_at_nodes() until no more small dangles are found */
 
997
         * and Vect_clean_small_angles_at_nodes() until no more small angles are found */
937
998
        do {
938
999
            G_message("%s", separator);
939
1000
            G_message(_("Break boundaries:"));
940
 
            Vect_break_lines(&Map, GV_BOUNDARY, NULL);
 
1001
            Vect_break_lines(&Tmp, GV_BOUNDARY, NULL);
941
1002
 
942
1003
            G_message("%s", separator);
943
1004
            G_message(_("Remove duplicates:"));
944
 
            Vect_remove_duplicates(&Map, GV_BOUNDARY, NULL);
 
1005
            Vect_remove_duplicates(&Tmp, GV_BOUNDARY, NULL);
945
1006
 
946
1007
            G_message("%s", separator);
947
1008
            G_message(_("Clean boundaries at nodes:"));
948
1009
            nmodif =
949
 
                Vect_clean_small_angles_at_nodes(&Map, GV_BOUNDARY, NULL);
 
1010
                Vect_clean_small_angles_at_nodes(&Tmp, GV_BOUNDARY, NULL);
950
1011
        } while (nmodif > 0);
951
1012
 
952
1013
        G_message("%s", separator);
953
1014
        if (type & GV_BOUNDARY) {       /* that means lines were converted boundaries */
954
1015
            G_message(_("Change boundary dangles to lines:"));
955
 
            Vect_chtype_dangles(&Map, -1.0, NULL);
 
1016
            Vect_chtype_dangles(&Tmp, -1.0, NULL);
956
1017
        }
957
1018
        else {
958
1019
            G_message(_("Change dangles to lines:"));
959
 
            Vect_remove_dangles(&Map, GV_BOUNDARY, -1.0, NULL);
 
1020
            Vect_remove_dangles(&Tmp, GV_BOUNDARY, -1.0, NULL);
960
1021
        }
961
1022
 
962
1023
        G_message("%s", separator);
963
1024
        if (type & GV_BOUNDARY) {
964
1025
            G_message(_("Change boundary bridges to lines:"));
965
 
            Vect_chtype_bridges(&Map, NULL);
 
1026
            Vect_chtype_bridges(&Tmp, NULL);
966
1027
        }
967
1028
        else {
968
1029
            G_message(_("Remove bridges:"));
969
 
            Vect_remove_bridges(&Map, NULL);
 
1030
            Vect_remove_bridges(&Tmp, NULL);
970
1031
        }
971
1032
 
 
1033
        /* merge boundaries */
 
1034
        G_message("%s", separator);
 
1035
        G_message(_("Merge boundaries:"));
 
1036
        Vect_merge_lines(&Tmp, GV_BOUNDARY, NULL, NULL);
 
1037
 
972
1038
        /* Boundaries are hopefully clean, build areas */
973
1039
        G_message("%s", separator);
974
 
        Vect_build_partial(&Map, GV_BUILD_ATTACH_ISLES);
 
1040
        Vect_build_partial(&Tmp, GV_BUILD_ATTACH_ISLES);
975
1041
 
976
1042
        /* Calculate new centroids for all areas, centroids have the same id as area */
977
 
        ncentr = Vect_get_num_areas(&Map);
 
1043
        ncentr = Vect_get_num_areas(&Tmp);
978
1044
        G_debug(3, "%d centroids/areas", ncentr);
979
1045
 
980
1046
        Centr = (CENTR *) G_calloc(ncentr + 1, sizeof(CENTR));
982
1048
        for (centr = 1; centr <= ncentr; centr++) {
983
1049
            Centr[centr].valid = 0;
984
1050
            Centr[centr].cats = Vect_new_cats_struct();
985
 
            ret = Vect_get_point_in_area(&Map, centr, &x, &y);
 
1051
            ret = Vect_get_point_in_area(&Tmp, centr, &x, &y);
986
1052
            if (ret < 0) {
987
1053
                G_warning(_("Cannot calculate area centroid"));
988
1054
                continue;
999
1065
 
1000
1066
        /* Go through all layers and find centroids for each polygon */
1001
1067
        for (layer = 0; layer < nlayers; layer++) {
1002
 
            G_message(_("Layer: %s"), layer_names[layer]);
 
1068
            G_message("%s", separator);
 
1069
            G_message(_("Find centroids for layer: %s"), layer_names[layer]);
1003
1070
            layer_id = layers[layer];
1004
1071
            Ogr_layer = OGR_DS_GetLayer(Ogr_ds, layer_id);
1005
1072
            OGR_L_ResetReading(Ogr_layer);
1019
1086
        }
1020
1087
 
1021
1088
        /* Write centroids */
 
1089
        G_message("%s", separator);
 
1090
        G_message(_("Write centroids:"));
 
1091
 
1022
1092
        n_overlaps = n_nocat = 0;
1023
1093
        total_area = overlap_area = nocat_area = 0.0;
1024
1094
        for (centr = 1; centr <= ncentr; centr++) {
1025
1095
            double area;
 
1096
            
 
1097
            G_percent(centr, ncentr, 2);
1026
1098
 
1027
 
            area = Vect_get_area_area(&Map, centr);
 
1099
            area = Vect_get_area_area(&Tmp, centr);
1028
1100
            total_area += area;
1029
1101
 
1030
1102
            if (!(Centr[centr].valid)) {
1050
1122
                otype = GV_POINT;
1051
1123
            else
1052
1124
                otype = GV_CENTROID;
1053
 
            Vect_write_line(&Map, otype, Points, Centr[centr].cats);
 
1125
            Vect_write_line(&Tmp, otype, Points, Centr[centr].cats);
1054
1126
        }
1055
1127
        if (Centr)
1056
1128
            G_free(Centr);
1057
 
        G_message("%s", separator);
1058
 
        Vect_build_partial(&Map, GV_BUILD_NONE);
1059
 
 
1060
 
        G_message("%s", separator);
1061
 
        Vect_build(&Map);
1062
 
 
1063
 
        G_message("%s", separator);
 
1129
            
 
1130
        Vect_spatial_index_destroy(&si);
 
1131
 
 
1132
        /* G_message("%s", separator); */
 
1133
        /* Vect_build_partial(&Map, GV_BUILD_NONE); */
 
1134
 
 
1135
        /* not needed */
 
1136
        /* G_message("%s", separator);
 
1137
           Vect_build(&Map); */
1064
1138
 
1065
1139
        if (n_overlaps > 0) {
1066
1140
            G_warning(_("%d areas represent more (overlapping) features, because polygons overlap "
1069
1143
                      n_overlaps, nlayers + 1);
1070
1144
        }
1071
1145
 
1072
 
        sprintf(buf, _("%d input polygons"), n_polygons);
1073
 
        G_message(buf);
1074
 
        Vect_hist_write(&Map, buf);
1075
 
 
1076
 
        sprintf(buf, _("Total area: %e (%d areas)"), total_area, ncentr);
1077
 
        G_message(buf);
1078
 
        Vect_hist_write(&Map, buf);
1079
 
 
1080
 
        sprintf(buf, _("Overlapping area: %e (%d areas)"), overlap_area,
 
1146
        G_message("%s", separator);
 
1147
 
 
1148
        Vect_hist_write(&Map, separator);
 
1149
        Vect_hist_write(&Map, "\n");
 
1150
        sprintf(buf, _("%d input polygons\n"), n_polygons);
 
1151
        G_message(_("%d input polygons"), n_polygons);
 
1152
        Vect_hist_write(&Map, buf);
 
1153
 
 
1154
        sprintf(buf, _("Total area: %G (%d areas)\n"), total_area, ncentr);
 
1155
        G_message(_("Total area: %G (%d areas)"), total_area, ncentr);
 
1156
        Vect_hist_write(&Map, buf);
 
1157
 
 
1158
        sprintf(buf, _("Overlapping area: %G (%d areas)\n"), overlap_area,
1081
1159
                n_overlaps);
1082
 
        G_message(buf);
 
1160
        G_message(_("Overlapping area: %G (%d areas)"), overlap_area,
 
1161
                  n_overlaps);
1083
1162
        Vect_hist_write(&Map, buf);
1084
1163
 
1085
 
        sprintf(buf, _("Area without category: %e (%d areas)"), nocat_area,
 
1164
        sprintf(buf, _("Area without category: %G (%d areas)\n"), nocat_area,
1086
1165
                n_nocat);
1087
 
        G_message(buf);
 
1166
        G_message(_("Area without category: %G (%d areas)"), nocat_area,
 
1167
                  n_nocat);
1088
1168
        Vect_hist_write(&Map, buf);
 
1169
        G_message("%s", separator);
1089
1170
    }
1090
1171
 
1091
1172
    /* needed?
1092
1173
     * OGR_DS_Destroy( Ogr_ds );
1093
1174
     */
1094
1175
 
 
1176
    /* Copy temporary vector to output vector */
 
1177
    Vect_copy_map_lines(&Tmp, &Map);
 
1178
    Vect_close(&Tmp);
 
1179
    Vect_delete(tempvect);
 
1180
 
 
1181
    Vect_build(&Map);
1095
1182
    Vect_close(&Map);
1096
1183
 
1097
1184