GRASS Programmer's Manual  6.4.2(2012)
vinput2d.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines