GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /*- 00003 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Fall 1993 00004 * University of Illinois 00005 * US Army Construction Engineering Research Lab 00006 * Copyright 1993, H. Mitasova (University of Illinois), 00007 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL) 00008 * 00009 * modified by McCauley in August 1995 00010 * modified by Mitasova in August 1995 00011 * modofied by Mitasova in Nov 1999 (dmax fix) 00012 * 00013 */ 00014 00015 #include <stdio.h> 00016 #include <stdlib.h> 00017 #include <math.h> 00018 #include <grass/bitmap.h> 00019 #include <grass/linkm.h> 00020 #include <grass/gis.h> 00021 #include <grass/dbmi.h> 00022 #include <grass/Vect.h> 00023 #include <grass/glocale.h> 00024 00025 #include <grass/interpf.h> 00026 00027 int IL_vector_input_data_2d(struct interp_params *params, struct Map_info *Map, /* input vector map */ 00028 /* as z values may be used: 1) z coordinates in 3D file -> field = 0 00029 * 2) categories -> field > 0, zcol = NULL 00030 * 3) attributes -> field > 0, zcol != NULL */ 00031 int field, /* category field number */ 00032 char *zcol, /* name of the column containing z values */ 00033 char *scol, /* name of the column containing smooth values */ 00034 struct tree_info *info, /* quadtree info */ 00035 double *xmin, double *xmax, double *ymin, double *ymax, double *zmin, double *zmax, int *n_points, /* number of points used for interpolation */ 00036 double *dmax) 00037 00038 /* 00039 * Inserts input data inside the region into a quad tree. Also translates 00040 * data. Returns number of segments in the quad tree. 00041 */ 00042 { 00043 double dmax2; /* max distance between points squared */ 00044 double c1, c2, c3, c4; 00045 int i, line, k = 0, nnodes; 00046 double ns_res, ew_res; 00047 int npoint, OUTRANGE; 00048 int totsegm; 00049 struct quaddata *data = (struct quaddata *)info->root->data; 00050 double xprev, yprev, zprev, x1, y1, z1, d1, xt, yt, z, sm; 00051 struct line_pnts *Points; 00052 struct line_cats *Cats; 00053 int times, j1, k1, ltype, cat, zctype = 0, sctype = 0; 00054 struct field_info *Fi; 00055 dbDriver *driver; 00056 dbHandle handle; 00057 dbString stmt; 00058 dbCatValArray zarray, sarray; 00059 00060 OUTRANGE = 0; 00061 npoint = 0; 00062 00063 G_debug(2, "IL_vector_input_data_2d(): field = %d, zcol = %s, scol = %s", 00064 field, zcol, scol); 00065 ns_res = (data->ymax - data->y_orig) / data->n_rows; 00066 ew_res = (data->xmax - data->x_orig) / data->n_cols; 00067 dmax2 = *dmax * *dmax; 00068 00069 Points = Vect_new_line_struct(); /* init line_pnts struct */ 00070 Cats = Vect_new_cats_struct(); 00071 00072 if (field == 0 && !Vect_is_3d(Map)) 00073 G_fatal_error(_("Vector map <%s> is not 3D"), Vect_get_full_name(Map)); 00074 00075 if (field > 0 && zcol != NULL) { /* open db driver */ 00076 G_verbose_message(_("Loading data from attribute table ...")); 00077 Fi = Vect_get_field(Map, field); 00078 if (Fi == NULL) 00079 G_fatal_error(_("Database connection not defined for layer %d"), 00080 field); 00081 G_debug(3, " driver = %s database = %s table = %s", Fi->driver, 00082 Fi->database, Fi->table); 00083 db_init_handle(&handle); 00084 db_init_string(&stmt); 00085 driver = db_start_driver(Fi->driver); 00086 db_set_handle(&handle, Fi->database, NULL); 00087 if (db_open_database(driver, &handle) != DB_OK) 00088 G_fatal_error(_("Unable to open database <%s> by driver <%s>"), 00089 Fi->database, Fi->driver); 00090 00091 zctype = db_column_Ctype(driver, Fi->table, zcol); 00092 G_debug(3, " zcol C type = %d", zctype); 00093 if (zctype == -1) 00094 G_fatal_error(_("Column <%s> not found"), 00095 zcol); 00096 if (zctype != DB_C_TYPE_INT && zctype != DB_C_TYPE_DOUBLE) 00097 G_fatal_error(_("Data type of column <%s> must be numeric"), zcol); 00098 00099 db_CatValArray_init(&zarray); 00100 G_debug(3, "RST SQL WHERE: %s", params->wheresql); 00101 db_select_CatValArray(driver, Fi->table, Fi->key, zcol, 00102 params->wheresql, &zarray); 00103 00104 if (scol != NULL) { 00105 sctype = db_column_Ctype(driver, Fi->table, scol); 00106 G_debug(3, " scol C type = %d", sctype); 00107 if (sctype == -1) 00108 G_fatal_error(_("Column <%s> not found"), scol); 00109 if (sctype != DB_C_TYPE_INT && sctype != DB_C_TYPE_DOUBLE) 00110 G_fatal_error(_("Data type of column <%s> must be numeric"), scol); 00111 00112 db_CatValArray_init(&sarray); 00113 db_select_CatValArray(driver, Fi->table, Fi->key, scol, 00114 params->wheresql, &sarray); 00115 } 00116 00117 db_close_database_shutdown_driver(driver); 00118 } 00119 00120 /* Lines without nodes */ 00121 G_message(_("Reading features from vector map ...")); 00122 sm = 0; 00123 line = 1; 00124 while ((ltype = Vect_read_next_line(Map, Points, Cats)) != -2) { 00125 00126 if (!(ltype & (GV_POINT | GV_LINE | GV_BOUNDARY))) 00127 continue; 00128 00129 if (field > 0) { /* use cat or attribute */ 00130 Vect_cat_get(Cats, field, &cat); 00131 00132 /* line++; */ 00133 if (zcol == NULL) { /* use categories */ 00134 z = (double)cat; 00135 } 00136 else { /* read att from db */ 00137 int ret, intval; 00138 00139 if (zctype == DB_C_TYPE_INT) { 00140 ret = db_CatValArray_get_value_int(&zarray, cat, &intval); 00141 z = intval; 00142 } 00143 else { /* DB_C_TYPE_DOUBLE */ 00144 ret = db_CatValArray_get_value_double(&zarray, cat, &z); 00145 } 00146 00147 if (ret != DB_OK) { 00148 if (params->wheresql != NULL) 00149 /* G_message(_("Database record for cat %d not used due to SQL statement")); */ 00150 /* do nothing in this case to not confuse user. Or implement second cat list */ 00151 ; 00152 else 00153 G_warning(_("Database record for cat %d not found"), 00154 cat); 00155 continue; 00156 } 00157 00158 if (scol != NULL) { 00159 if (sctype == DB_C_TYPE_INT) { 00160 ret = 00161 db_CatValArray_get_value_int(&sarray, cat, 00162 &intval); 00163 sm = intval; 00164 } 00165 else { /* DB_C_TYPE_DOUBLE */ 00166 ret = 00167 db_CatValArray_get_value_double(&sarray, cat, 00168 &sm); 00169 } 00170 if (sm < 0.0) 00171 G_fatal_error(_("Negative value of smoothing detected: sm must be >= 0")); 00172 } 00173 G_debug(5, " z = %f sm = %f", z, sm); 00174 } 00175 } 00176 00177 if (Vect_level(Map) == 1) { 00178 /* Insert all points including nodes (end points) */ 00179 for (i = 0; i < Points->n_points; i++) { 00180 if (field == 0) 00181 z = Points->z[i]; 00182 process_point(Points->x[i], Points->y[i], z, sm, info, 00183 params->zmult, xmin, xmax, ymin, ymax, zmin, 00184 zmax, &npoint, &OUTRANGE, &k); 00185 00186 } 00187 } 00188 else { 00189 /* Insert all points except nodes (end points) */ 00190 for (i = 1; i < Points->n_points - 1; i++) { 00191 if (field == 0) 00192 z = Points->z[i]; 00193 process_point(Points->x[i], Points->y[i], z, sm, info, 00194 params->zmult, xmin, xmax, ymin, ymax, zmin, 00195 zmax, &npoint, &OUTRANGE, &k); 00196 00197 } 00198 } 00199 00200 /* Check all segments */ 00201 xprev = Points->x[0]; 00202 yprev = Points->y[0]; 00203 zprev = Points->z[0]; 00204 for (i = 1; i < Points->n_points; i++) { 00205 /* compare the distance between current and previous */ 00206 x1 = Points->x[i]; 00207 y1 = Points->y[i]; 00208 z1 = Points->z[i]; 00209 00210 xt = x1 - xprev; 00211 yt = y1 - yprev; 00212 d1 = (xt * xt + yt * yt); 00213 if ((d1 > dmax2) && (dmax2 != 0.)) { 00214 times = (int)(d1 / dmax2 + 0.5); 00215 for (j1 = 0; j1 < times; j1++) { 00216 xt = x1 - j1 * ((x1 - xprev) / times); 00217 yt = y1 - j1 * ((y1 - yprev) / times); 00218 if (field == 0) 00219 z = z1 - j1 * ((z1 - zprev) / times); 00220 00221 process_point(xt, yt, z, sm, info, params->zmult, 00222 xmin, xmax, ymin, ymax, zmin, zmax, &npoint, 00223 &OUTRANGE, &k); 00224 } 00225 } 00226 xprev = x1; 00227 yprev = y1; 00228 zprev = z1; 00229 } 00230 } 00231 00232 /* Process all nodes */ 00233 G_message(_("Reading nodes from vector map ...")); 00234 nnodes = Vect_get_num_nodes(Map); 00235 for (k1 = 1; k1 <= nnodes; k1++) { 00236 G_debug(5, " node %d", k1); 00237 G_percent(k1, nnodes, 1); 00238 Vect_get_node_coor(Map, k1, &x1, &y1, &z); 00239 00240 /* TODO: check more lines ? */ 00241 if (field > 0) { 00242 line = abs(Vect_get_node_line(Map, k1, 0)); 00243 ltype = Vect_read_line(Map, NULL, Cats, line); 00244 Vect_cat_get(Cats, field, &cat); 00245 if (zcol == NULL) { /* use categories */ 00246 if (cat < 0) 00247 continue; 00248 z = (double)cat; 00249 } 00250 else { /* read att from db */ 00251 int ret, intval; 00252 00253 if (cat == 0) 00254 continue; 00255 00256 if (zctype == DB_C_TYPE_INT) { 00257 ret = db_CatValArray_get_value_int(&zarray, cat, &intval); 00258 z = intval; 00259 } 00260 else { /* DB_C_TYPE_DOUBLE */ 00261 ret = db_CatValArray_get_value_double(&zarray, cat, &z); 00262 } 00263 00264 if (ret != DB_OK) { 00265 if (params->wheresql != NULL) 00266 /* G_message(_("Database record for cat %d not used due to SQL statement")); */ 00267 /* do nothing in this case to not confuse user. Or implement second cat list */ 00268 ; 00269 else 00270 G_warning(_("Database record for cat %d not found"), 00271 cat); 00272 continue; 00273 } 00274 00275 if (scol != NULL) { 00276 if (sctype == DB_C_TYPE_INT) { 00277 ret = 00278 db_CatValArray_get_value_int(&sarray, cat, 00279 &intval); 00280 sm = intval; 00281 } 00282 else { /* DB_C_TYPE_DOUBLE */ 00283 ret = 00284 db_CatValArray_get_value_double(&sarray, cat, 00285 &sm); 00286 } 00287 if (sm < 0.0) 00288 G_fatal_error(_("Negative value of smoothing detected: sm must be >= 0")); 00289 } 00290 G_debug(5, " z = %f sm = %f", z, sm); 00291 } 00292 } 00293 00294 process_point(x1, y1, z, sm, info, params->zmult, xmin, xmax, ymin, 00295 ymax, zmin, zmax, &npoint, &OUTRANGE, &k); 00296 } 00297 00298 if (field > 0 && zcol != NULL) 00299 db_CatValArray_free(&zarray); 00300 if (scol != NULL) { 00301 db_CatValArray_free(&sarray); 00302 } 00303 00304 c1 = *xmin - data->x_orig; 00305 c2 = data->xmax - *xmax; 00306 c3 = *ymin - data->y_orig; 00307 c4 = data->ymax - *ymax; 00308 if ((c1 > 5 * ew_res) || (c2 > 5 * ew_res) || (c3 > 5 * ns_res) || 00309 (c4 > 5 * ns_res)) { 00310 static int once = 0; 00311 00312 if (!once) { 00313 once = 1; 00314 G_warning(_("Strip exists with insufficient data")); 00315 } 00316 } 00317 00318 totsegm = 00319 translate_quad(info->root, data->x_orig, data->y_orig, *zmin, 4); 00320 if (!totsegm) 00321 return 0; 00322 data->x_orig = 0; 00323 data->y_orig = 0; 00324 00325 /* G_read_vector_timestamp(name,mapset,ts); */ 00326 00327 if (OUTRANGE > 0) 00328 G_warning(_("There are points outside specified 2D/3D region - %d points ignored"), 00329 OUTRANGE); 00330 if (npoint > 0) 00331 G_important_message(_("Ignoring %d points (too dense)"), npoint); 00332 npoint = k - npoint - OUTRANGE; 00333 if (npoint < params->kmin) { 00334 if (npoint != 0) { 00335 G_warning(_("%d points given for interpolation (after thinning) is less than given NPMIN=%d"), 00336 npoint, params->kmin); 00337 params->kmin = npoint; 00338 } 00339 else { 00340 G_warning(_("Zero points in the given region")); 00341 return -1; 00342 } 00343 } 00344 if (npoint > params->KMAX2 && params->kmin <= params->kmax) { 00345 G_warning(_("Segmentation parameters set to invalid values: npmin= %d, segmax= %d " 00346 "for smooth connection of segments, npmin > segmax (see manual)"), 00347 params->kmin, params->kmax); 00348 return -1; 00349 } 00350 if (npoint < params->KMAX2 && params->kmax != params->KMAX2) 00351 G_warning(_("There are less than %d points for interpolation. No " 00352 "segmentation is necessary, to run the program faster set " 00353 "segmax=%d (see manual)"), params->KMAX2, params->KMAX2); 00354 00355 G_message(_("Number of points from vector map %d"), k); 00356 G_verbose_message(_("Number of points outside of 2D/3D region %d"), 00357 OUTRANGE); 00358 G_message(_("Number of points being used %d"), npoint); 00359 00360 *n_points = npoint; 00361 return (totsegm); 00362 } 00363 00364 int process_point(double x, double y, double z, double sm, struct tree_info *info, /* quadtree info */ 00365 double zmult, /* multiplier for z-values */ 00366 double *xmin, 00367 double *xmax, 00368 double *ymin, 00369 double *ymax, 00370 double *zmin, 00371 double *zmax, int *npoint, int *OUTRANGE, int *total) 00372 { 00373 struct triple *point; 00374 double c1, c2, c3, c4; 00375 int a; 00376 static int first_time = 1; 00377 struct quaddata *data = (struct quaddata *)info->root->data; 00378 00379 00380 (*total)++; 00381 00382 00383 z = z * zmult; 00384 c1 = x - data->x_orig; 00385 c2 = data->xmax - x; 00386 c3 = y - data->y_orig; 00387 c4 = data->ymax - y; 00388 00389 if (!((c1 >= 0) && (c2 >= 0) && (c3 >= 0) && (c4 >= 0))) { 00390 if (!(*OUTRANGE)) { 00391 G_warning(_("Some points outside of region (ignored)")); 00392 } 00393 (*OUTRANGE)++; 00394 } 00395 else { 00396 if (!(point = quad_point_new(x, y, z, sm))) { 00397 G_warning(_("Unable to allocate memory")); 00398 return -1; 00399 } 00400 a = MT_insert(point, info, info->root, 4); 00401 if (a == 0) { 00402 (*npoint)++; 00403 } 00404 if (a < 0) { 00405 G_warning(_("Unable to insert %f,%f,%f a = %d"), x, y, z, a); 00406 return -1; 00407 } 00408 free(point); 00409 if (first_time) { 00410 first_time = 0; 00411 *xmin = x; 00412 *ymin = y; 00413 *zmin = z; 00414 *xmax = x; 00415 *ymax = y; 00416 *zmax = z; 00417 } 00418 *xmin = amin1(*xmin, x); 00419 *ymin = amin1(*ymin, y); 00420 *zmin = amin1(*zmin, z); 00421 *xmax = amax1(*xmax, x); 00422 *ymax = amax1(*ymax, y); 00423 *zmax = amax1(*zmax, z); 00424 } 00425 return 1; 00426 }