GRASS Programmer's Manual  6.4.2(2012)
N_arrays_io.c
Go to the documentation of this file.
00001 
00002 /*****************************************************************************
00003 *
00004 * MODULE:       Grass PDE Numerical Library
00005 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
00006 *               soerengebbert <at> gmx <dot> de
00007 *               
00008 * PURPOSE:      IO array managment functions 
00009 *               part of the gpde library
00010 *
00011 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
00012 *
00013 *               This program is free software under the GNU General Public
00014 *               License (>=v2). Read the file COPYING that comes with GRASS
00015 *               for details.
00016 *
00017 *****************************************************************************/
00018 
00019 #include <math.h>
00020 
00021 #include "grass/N_pde.h"
00022 #include "grass/glocale.h"
00023 
00024 /* ******************** 2D ARRAY FUNCTIONS *********************** */
00025 
00046 N_array_2d *N_read_rast_to_array_2d(char *name, N_array_2d * array)
00047 {
00048     int map;                    /*The rastermap */
00049     int x, y, cols, rows, type;
00050     void *rast;
00051     void *ptr;
00052     struct Cell_head region;
00053     N_array_2d *data = array;
00054 
00055     if (NULL == G_find_cell2(name, ""))
00056         G_fatal_error(_("Raster map <%s> not found"), name);
00057 
00058     /* Get the active region */
00059     G_get_set_window(&region);
00060 
00061     /*set the rows and cols */
00062     rows = region.rows;
00063     cols = region.cols;
00064 
00065     /*open the raster map */
00066     map = G_open_cell_old(name, G_find_cell2(name, ""));
00067     if (map < 0)
00068         G_fatal_error(_("Unable to open raster map <%s>"), name);
00069 
00070     type = G_get_raster_map_type(map);
00071 
00072     /*if the array is NULL create a new one with the data type of the raster map */
00073     /*the offset is 0 by default */
00074     if (data == NULL) {
00075         if (type == DCELL_TYPE) {
00076             data = N_alloc_array_2d(cols, rows, 0, DCELL_TYPE);
00077         }
00078         if (type == FCELL_TYPE) {
00079             data = N_alloc_array_2d(cols, rows, 0, FCELL_TYPE);
00080         }
00081         if (type == CELL_TYPE) {
00082             data = N_alloc_array_2d(cols, rows, 0, CELL_TYPE);
00083         }
00084     }
00085     else {
00086         /*Check the array sizes */
00087         if (data->cols != cols)
00088             G_fatal_error
00089                 ("N_read_rast_to_array_2d: the data array size is different from the current region settings");
00090         if (data->rows != rows)
00091             G_fatal_error
00092                 ("N_read_rast_to_array_2d: the data array size is different from the current region settings");
00093     }
00094 
00095     rast = G_allocate_raster_buf(type);
00096 
00097     G_message(_("Reading raster map <%s> into memory"), name);
00098 
00099     for (y = 0; y < rows; y++) {
00100         G_percent(y, rows - 1, 10);
00101 
00102         if (!G_get_raster_row(map, rast, y, type)) {
00103             G_close_cell(map);
00104             G_fatal_error(_("Could not get raster row"));
00105         }
00106 
00107         for (x = 0, ptr = rast; x < cols;
00108              x++, ptr = G_incr_void_ptr(ptr, G_raster_size(type))) {
00109             if (type == CELL_TYPE) {
00110                 if (G_is_c_null_value(ptr)) {
00111                     N_put_array_2d_value_null(data, x, y);
00112                 }
00113                 else {
00114                     if (data->type == CELL_TYPE)
00115                         N_put_array_2d_c_value(data, x, y,
00116                                                (CELL) * (CELL *) ptr);
00117                     if (data->type == FCELL_TYPE)
00118                         N_put_array_2d_f_value(data, x, y,
00119                                                (FCELL) * (CELL *) ptr);
00120                     if (data->type == DCELL_TYPE)
00121                         N_put_array_2d_d_value(data, x, y,
00122                                                (DCELL) * (CELL *) ptr);
00123                 }
00124             }
00125             if (type == FCELL_TYPE) {
00126                 if (G_is_f_null_value(ptr)) {
00127                     N_put_array_2d_value_null(data, x, y);
00128                 }
00129                 else {
00130                     if (data->type == CELL_TYPE)
00131                         N_put_array_2d_c_value(data, x, y,
00132                                                (CELL) * (FCELL *) ptr);
00133                     if (data->type == FCELL_TYPE)
00134                         N_put_array_2d_f_value(data, x, y,
00135                                                (FCELL) * (FCELL *) ptr);
00136                     if (data->type == DCELL_TYPE)
00137                         N_put_array_2d_d_value(data, x, y,
00138                                                (DCELL) * (FCELL *) ptr);
00139                 }
00140             }
00141             if (type == DCELL_TYPE) {
00142                 if (G_is_d_null_value(ptr)) {
00143                     N_put_array_2d_value_null(data, x, y);
00144                 }
00145                 else {
00146                     if (data->type == CELL_TYPE)
00147                         N_put_array_2d_c_value(data, x, y,
00148                                                (CELL) * (DCELL *) ptr);
00149                     if (data->type == FCELL_TYPE)
00150                         N_put_array_2d_f_value(data, x, y,
00151                                                (FCELL) * (DCELL *) ptr);
00152                     if (data->type == DCELL_TYPE)
00153                         N_put_array_2d_d_value(data, x, y,
00154                                                (DCELL) * (DCELL *) ptr);
00155                 }
00156             }
00157         }
00158     }
00159 
00160     /* Close file */
00161     if (G_close_cell(map) < 0)
00162         G_fatal_error(_("Unable to close input map"));
00163 
00164     return data;
00165 }
00166 
00181 void N_write_array_2d_to_rast(N_array_2d * array, char *name)
00182 {
00183     int map;                    /*The rastermap */
00184     int x, y, cols, rows, count, type;
00185     CELL *rast = NULL;
00186     FCELL *frast = NULL;
00187     DCELL *drast = NULL;
00188     struct Cell_head region;
00189 
00190     if (!array)
00191         G_fatal_error(_("N_array_2d * array is empty"));
00192 
00193     /* Get the current region */
00194     G_get_set_window(&region);
00195 
00196     rows = region.rows;
00197     cols = region.cols;
00198     type = array->type;
00199 
00200     /*Open the new map */
00201     map = G_open_raster_new(name, type);
00202     if (map < 0)
00203         G_fatal_error(_("Unable to create raster map <%s>"), name);
00204 
00205     if (type == CELL_TYPE)
00206         rast = G_allocate_raster_buf(type);
00207     if (type == FCELL_TYPE)
00208         frast = G_allocate_raster_buf(type);
00209     if (type == DCELL_TYPE)
00210         drast = G_allocate_raster_buf(type);
00211 
00212     G_message(_("Write 2d array to raster map <%s>"), name);
00213 
00214     count = 0;
00215     for (y = 0; y < rows; y++) {
00216         G_percent(y, rows - 1, 10);
00217         for (x = 0; x < cols; x++) {
00218             if (type == CELL_TYPE)
00219                 rast[x] = N_get_array_2d_c_value(array, x, y);
00220             if (type == FCELL_TYPE)
00221                 frast[x] = N_get_array_2d_f_value(array, x, y);
00222             if (type == DCELL_TYPE)
00223                 drast[x] = N_get_array_2d_d_value(array, x, y);
00224         }
00225         if (type == CELL_TYPE)
00226             if (!G_put_c_raster_row(map, rast)) {
00227                 G_unopen_cell(map);     /*unopen the new raster map */
00228                 G_fatal_error(_("Unable to write raster row %i"), y);
00229             }
00230         if (type == FCELL_TYPE)
00231             if (!G_put_f_raster_row(map, frast)) {
00232                 G_unopen_cell(map);     /*unopen the new raster map */
00233                 G_fatal_error(_("Unable to write raster row %i"), y);
00234             }
00235         if (type == DCELL_TYPE)
00236             if (!G_put_d_raster_row(map, drast)) {
00237                 G_unopen_cell(map);     /*unopen the new raster map */
00238                 G_fatal_error(_("Unable to write raster row %i"), y);
00239             }
00240     }
00241 
00242     /* Close file */
00243     if (G_close_cell(map) < 0)
00244         G_fatal_error(_("Unable to close input map"));
00245 
00246     return;
00247 
00248 }
00249 
00250 
00251 /* ******************** 3D ARRAY FUNCTIONS *********************** */
00252 
00276 N_array_3d *N_read_rast3d_to_array_3d(char *name, N_array_3d * array,
00277                                       int mask)
00278 {
00279     void *map = NULL;           /*The 3D Rastermap */
00280     int changemask = 0;
00281     int x, y, z, cols, rows, depths, type;
00282     double d1 = 0, f1 = 0;
00283     N_array_3d *data = array;
00284     G3D_Region region;
00285 
00286 
00287     /*get the current region */
00288     G3d_getWindow(&region);
00289 
00290     cols = region.cols;
00291     rows = region.rows;
00292     depths = region.depths;
00293 
00294 
00295     if (NULL == G_find_grid3(name, ""))
00296         G3d_fatalError(_("3D raster map <%s> not found"), name);
00297 
00298     /*Open all maps with default region */
00299     map =
00300         G3d_openCellOld(name, G_find_grid3(name, ""), G3D_DEFAULT_WINDOW,
00301                         G3D_TILE_SAME_AS_FILE, G3D_USE_CACHE_DEFAULT);
00302 
00303     if (map == NULL)
00304         G3d_fatalError(_("Unable to open 3D raster map <%s>"), name);
00305 
00306     type = G3d_tileTypeMap(map);
00307 
00308     /*if the array is NULL create a new one with the data type of the volume map */
00309     /*the offset is 0 by default */
00310     if (data == NULL) {
00311         if (type == FCELL_TYPE) {
00312             data = N_alloc_array_3d(cols, rows, depths, 0, FCELL_TYPE);
00313         }
00314         if (type == DCELL_TYPE) {
00315             data = N_alloc_array_3d(cols, rows, depths, 0, DCELL_TYPE);
00316         }
00317     }
00318     else {
00319         /*Check the array sizes */
00320         if (data->cols != cols)
00321             G_fatal_error
00322                 ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
00323         if (data->rows != rows)
00324             G_fatal_error
00325                 ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
00326         if (data->depths != depths)
00327             G_fatal_error
00328                 ("N_read_rast_to_array_3d: the data array size is different from the current region settings");
00329     }
00330 
00331 
00332     G_message(_("Read g3d map <%s> into the memory"), name);
00333 
00334     /*if requested set the Mask on */
00335     if (mask) {
00336         if (G3d_maskFileExists()) {
00337             changemask = 0;
00338             if (G3d_maskIsOff(map)) {
00339                 G3d_maskOn(map);
00340                 changemask = 1;
00341             }
00342         }
00343     }
00344 
00345     for (z = 0; z < depths; z++) {      /*From the bottom to the top */
00346         G_percent(z, depths - 1, 10);
00347         for (y = 0; y < rows; y++) {
00348             for (x = 0; x < cols; x++) {
00349                 if (type == FCELL_TYPE) {
00350                     G3d_getValue(map, x, y, z, &f1, type);
00351                     if (G_is_f_null_value((void *)&f1)) {
00352                         N_put_array_3d_value_null(data, x, y, z);
00353                     }
00354                     else {
00355                         if (data->type == FCELL_TYPE)
00356                             N_put_array_3d_f_value(data, x, y, z, f1);
00357                         if (data->type == DCELL_TYPE)
00358                             N_put_array_3d_d_value(data, x, y, z, (double)f1);
00359                     }
00360                 }
00361                 else {
00362                     G3d_getValue(map, x, y, z, &d1, type);
00363                     if (G_is_d_null_value((void *)&d1)) {
00364                         N_put_array_3d_value_null(data, x, y, z);
00365                     }
00366                     else {
00367                         if (data->type == FCELL_TYPE)
00368                             N_put_array_3d_f_value(data, x, y, z, (float)d1);
00369                         if (data->type == DCELL_TYPE)
00370                             N_put_array_3d_d_value(data, x, y, z, d1);
00371                     }
00372 
00373                 }
00374             }
00375         }
00376     }
00377 
00378     /*We set the Mask off, if it was off before */
00379     if (mask) {
00380         if (G3d_maskFileExists())
00381             if (G3d_maskIsOn(map) && changemask)
00382                 G3d_maskOff(map);
00383     }
00384 
00385     /* Close files and exit */
00386     if (!G3d_closeCell(map))
00387         G3d_fatalError(map, NULL, 0, _("Error closing g3d file"));
00388 
00389     return data;
00390 }
00391 
00408 void N_write_array_3d_to_rast3d(N_array_3d * array, char *name, int mask)
00409 {
00410     void *map = NULL;           /*The 3D Rastermap */
00411     int changemask = 0;
00412     int x, y, z, cols, rows, depths, count, type;
00413     double d1 = 0.0, f1 = 0.0;
00414     N_array_3d *data = array;
00415     G3D_Region region;
00416 
00417     /*get the current region */
00418     G3d_getWindow(&region);
00419 
00420 
00421     cols = region.cols;
00422     rows = region.rows;
00423     depths = region.depths;
00424     type = data->type;
00425 
00426     /*Check the array sizes */
00427     if (data->cols != cols)
00428         G_fatal_error
00429             ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
00430     if (data->rows != rows)
00431         G_fatal_error
00432             ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
00433     if (data->depths != depths)
00434         G_fatal_error
00435             ("N_write_array_3d_to_rast3d: the data array size is different from the current region settings");
00436 
00437     /*Open the new map */
00438     if (type == DCELL_TYPE)
00439         map =
00440             G3d_openCellNew(name, DCELL_TYPE, G3D_USE_CACHE_DEFAULT, &region);
00441     else if (type == FCELL_TYPE)
00442         map =
00443             G3d_openCellNew(name, FCELL_TYPE, G3D_USE_CACHE_DEFAULT, &region);
00444 
00445     if (map == NULL)
00446         G3d_fatalError(_("Error opening g3d map <%s>"), name);
00447 
00448     G_message(_("Write 3d array to g3d map <%s>"), name);
00449 
00450     /*if requested set the Mask on */
00451     if (mask) {
00452         if (G3d_maskFileExists()) {
00453             changemask = 0;
00454             if (G3d_maskIsOff(map)) {
00455                 G3d_maskOn(map);
00456                 changemask = 1;
00457             }
00458         }
00459     }
00460 
00461     count = 0;
00462     for (z = 0; z < depths; z++) {      /*From the bottom to the top */
00463         G_percent(z, depths - 1, 10);
00464         for (y = 0; y < rows; y++) {
00465             for (x = 0; x < cols; x++) {
00466                 if (type == FCELL_TYPE) {
00467                     f1 = N_get_array_3d_f_value(data, x, y, z);
00468                     G3d_putFloat(map, x, y, z, f1);
00469                 }
00470                 else if (type == DCELL_TYPE) {
00471                     d1 = N_get_array_3d_d_value(data, x, y, z);
00472                     G3d_putDouble(map, x, y, z, d1);
00473                 }
00474             }
00475         }
00476     }
00477 
00478     /*We set the Mask off, if it was off before */
00479     if (mask) {
00480         if (G3d_maskFileExists())
00481             if (G3d_maskIsOn(map) && changemask)
00482                 G3d_maskOff(map);
00483     }
00484 
00485     /* Close files and exit */
00486     if (!G3d_closeCell(map))
00487         G3d_fatalError(map, NULL, 0, _("Error closing g3d file"));
00488 
00489     return;
00490 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines