GRASS Programmer's Manual
6.4.2(2012)
|
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(®ion); 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(®ion); 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(®ion); 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(®ion); 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, ®ion); 00441 else if (type == FCELL_TYPE) 00442 map = 00443 G3d_openCellNew(name, FCELL_TYPE, G3D_USE_CACHE_DEFAULT, ®ion); 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 }