GRASS Programmer's Manual
6.4.2(2012)
|
00001 00019 #include <stdlib.h> 00020 00021 #include <grass/gis.h> 00022 #include <grass/Vect.h> 00023 #include <grass/glocale.h> 00024 #include <grass/gstypes.h> 00025 00026 /* 00027 #define TRAK_MEM 00028 */ 00029 00030 #ifdef TRAK_MEM 00031 static int Tot_mem = 0; 00032 #endif 00033 00045 geoline *Gv_load_vect(const char *grassname, int *nlines) 00046 { 00047 struct Map_info map; 00048 struct line_pnts *points; 00049 geoline *top, *gln, *prev; 00050 int np, i, n, nareas, nl = 0, area, type, is3d; 00051 struct Cell_head wind; 00052 float vect[2][3]; 00053 const char *mapset; 00054 00055 mapset = G_find_vector2(grassname, ""); 00056 if (!mapset) { 00057 G_warning(_("Vector map <%s> not found"), grassname); 00058 return NULL; 00059 } 00060 00061 Vect_set_open_level(2); 00062 if (Vect_open_old(&map, grassname, "") == -1) { 00063 G_warning(_("Unable to open vector map <%s>"), 00064 G_fully_qualified_name(grassname, mapset)); 00065 return NULL; 00066 } 00067 00068 top = gln = (geoline *) G_malloc(sizeof(geoline)); /* G_fatal_error */ 00069 if (!top) { 00070 return NULL; 00071 } 00072 00073 prev = top; 00074 00075 #ifdef TRAK_MEM 00076 Tot_mem += sizeof(geoline); 00077 #endif 00078 00079 points = Vect_new_line_struct(); 00080 00081 G_get_set_window(&wind); 00082 Vect_set_constraint_region(&map, wind.north, wind.south, wind.east, 00083 wind.west, PORT_DOUBLE_MAX, -PORT_DOUBLE_MAX); 00084 00085 is3d = Vect_is_3d(&map); 00086 00087 /* Read areas */ 00088 n = Vect_get_num_areas(&map); 00089 nareas = 0; 00090 G_debug(3, "Reading vector areas (nareas = %d)", n); 00091 for (area = 1; area <= n; area++) { 00092 G_debug(3, " area %d", area); 00093 00094 Vect_get_area_points(&map, area, points); 00095 if (points->n_points < 3) 00096 continue; 00097 gln->type = OGSF_POLYGON; 00098 gln->npts = np = points->n_points; 00099 G_debug(3, " np = %d", np); 00100 00101 if (is3d) { 00102 gln->dims = 3; 00103 gln->p3 = (Point3 *) G_calloc(np, sizeof(Point3)); /* G_fatal_error */ 00104 if (!gln->p3) { 00105 return (NULL); 00106 } 00107 #ifdef TRAK_MEM 00108 Tot_mem += (np * sizeof(Point3)); 00109 #endif 00110 } 00111 else { 00112 gln->dims = 2; 00113 gln->p2 = (Point2 *) G_calloc(np, sizeof(Point2)); /* G_fatal_error */ 00114 if (!gln->p2) { 00115 return (NULL); 00116 } 00117 #ifdef TRAK_MEM 00118 Tot_mem += (np * sizeof(Point2)); 00119 #endif 00120 } 00121 00122 for (i = 0; i < np; i++) { 00123 if (is3d) { 00124 gln->p3[i][X] = points->x[i]; 00125 gln->p3[i][Y] = points->y[i]; 00126 gln->p3[i][Z] = points->z[i]; 00127 } 00128 else { 00129 gln->p2[i][X] = points->x[i]; 00130 gln->p2[i][Y] = points->y[i]; 00131 } 00132 } 00133 /* Calc normal (should be average) */ 00134 if (is3d) { 00135 vect[0][X] = (float)(gln->p3[0][X] - gln->p3[1][X]); 00136 vect[0][Y] = (float)(gln->p3[0][Y] - gln->p3[1][Y]); 00137 vect[0][Z] = (float)(gln->p3[0][Z] - gln->p3[1][Z]); 00138 vect[1][X] = (float)(gln->p3[2][X] - gln->p3[1][X]); 00139 vect[1][Y] = (float)(gln->p3[2][Y] - gln->p3[1][Y]); 00140 vect[1][Z] = (float)(gln->p3[2][Z] - gln->p3[1][Z]); 00141 GS_v3cross(vect[1], vect[0], gln->norm); 00142 00143 } 00144 00145 gln->next = (geoline *) G_malloc(sizeof(geoline)); /* G_fatal_error */ 00146 if (!gln->next) { 00147 return (NULL); 00148 } 00149 00150 #ifdef TRAK_MEM 00151 Tot_mem += sizeof(geoline); 00152 #endif 00153 00154 prev = gln; 00155 gln = gln->next; 00156 nareas++; 00157 } 00158 G_debug(3, "%d areas loaded", nareas); 00159 00160 /* Read all lines */ 00161 G_debug(3, "Reading vector lines ..."); 00162 while (-1 < (type = Vect_read_next_line(&map, points, NULL))) { 00163 G_debug(3, "line type = %d", type); 00164 00165 if (type & (GV_LINES | GV_FACE)) { 00166 if (type & (GV_LINES)) { 00167 gln->type = OGSF_LINE; 00168 } 00169 else { 00170 gln->type = OGSF_POLYGON; 00171 /* Vect_append_point ( points, points->x[0], points->y[0], points->z[0] ); */ 00172 } 00173 00174 gln->npts = np = points->n_points; 00175 G_debug(3, " np = %d", np); 00176 00177 if (is3d) { 00178 gln->dims = 3; 00179 gln->p3 = (Point3 *) G_calloc(np, sizeof(Point3)); /* G_fatal_error */ 00180 if (!gln->p3) { 00181 return (NULL); 00182 } 00183 #ifdef TRAK_MEM 00184 Tot_mem += (np * sizeof(Point3)); 00185 #endif 00186 } 00187 else { 00188 gln->dims = 2; 00189 gln->p2 = (Point2 *) G_calloc(np, sizeof(Point2)); /* G_fatal_error */ 00190 if (!gln->p2) { 00191 return (NULL); 00192 } 00193 #ifdef TRAK_MEM 00194 Tot_mem += (np * sizeof(Point2)); 00195 #endif 00196 } 00197 00198 for (i = 0; i < np; i++) { 00199 if (is3d) { 00200 gln->p3[i][X] = points->x[i]; 00201 gln->p3[i][Y] = points->y[i]; 00202 gln->p3[i][Z] = points->z[i]; 00203 } 00204 else { 00205 gln->p2[i][X] = points->x[i]; 00206 gln->p2[i][Y] = points->y[i]; 00207 } 00208 } 00209 /* Calc normal (should be average) */ 00210 if (is3d && gln->type == OGSF_POLYGON) { 00211 vect[0][X] = (float)(gln->p3[0][X] - gln->p3[1][X]); 00212 vect[0][Y] = (float)(gln->p3[0][Y] - gln->p3[1][Y]); 00213 vect[0][Z] = (float)(gln->p3[0][Z] - gln->p3[1][Z]); 00214 vect[1][X] = (float)(gln->p3[2][X] - gln->p3[1][X]); 00215 vect[1][Y] = (float)(gln->p3[2][Y] - gln->p3[1][Y]); 00216 vect[1][Z] = (float)(gln->p3[2][Z] - gln->p3[1][Z]); 00217 GS_v3cross(vect[1], vect[0], gln->norm); 00218 G_debug(3, "norm %f %f %f", gln->norm[0], gln->norm[1], 00219 gln->norm[2]); 00220 } 00221 00222 gln->next = (geoline *) G_malloc(sizeof(geoline)); /* G_fatal_error */ 00223 if (!gln->next) { 00224 return (NULL); 00225 } 00226 #ifdef TRAK_MEM 00227 Tot_mem += sizeof(geoline); 00228 #endif 00229 00230 prev = gln; 00231 gln = gln->next; 00232 nl++; 00233 } 00234 } 00235 G_debug(3, "%d lines loaded", nl); 00236 00237 nl += nareas; 00238 00239 prev->next = NULL; 00240 G_free(gln); 00241 00242 #ifdef TRAK_MEM 00243 Tot_mem -= sizeof(geoline); 00244 #endif 00245 00246 Vect_close(&map); 00247 00248 if (!nl) { 00249 G_warning(_("No features from vector map <%s> fall within current region"), 00250 G_fully_qualified_name(grassname, mapset)); 00251 return (NULL); 00252 } 00253 else { 00254 G_message(_("Vector map <%s> loaded (%d features)"), 00255 G_fully_qualified_name(grassname, mapset), nl); 00256 } 00257 00258 *nlines = nl; 00259 00260 #ifdef TRAK_MEM 00261 G_debug(3, "Total vect memory = %d Kbytes", Tot_mem / 1000); 00262 #endif 00263 00264 return (top); 00265 } 00266 00272 void sub_Vectmem(int minus) 00273 { 00274 #ifdef TRAK_MEM 00275 { 00276 Tot_mem -= minus; 00277 } 00278 #endif 00279 00280 return; 00281 }