openModeller
Version 1.4.0
|
00001 00030 #include <openmodeller/env_io/GdalRaster.hh> 00031 #include <openmodeller/env_io/Raster.hh> 00032 #include <openmodeller/env_io/Map.hh> 00033 #include <openmodeller/env_io/GeoTransform.hh> 00034 #include <openmodeller/Log.hh> 00035 #include <openmodeller/Exceptions.hh> 00036 #include <openmodeller/MapFormat.hh> 00037 00038 #include <gdal_version.h> 00039 #include <gdal.h> 00040 #include <gdal_priv.h> 00041 #include <cpl_string.h> 00042 #include <gdalwarper.h> 00043 00044 #include <string.h> 00045 00046 using std::string; 00047 using std::vector; 00048 00049 #include <utility> 00050 using std::pair; 00051 00052 #ifdef MPI_FOUND 00053 #include "mpi.h" 00054 #endif 00055 00056 struct GDAL_Format 00057 { 00058 char const *GDalDriverName; 00059 GDALDataType dataType; 00060 bool hasMeta; 00061 }; 00062 00063 GDAL_Format Formats[8] = 00064 { 00065 // Floating GeoTiff 00066 { "GTiff", 00067 (sizeof(Scalar) == 4) ? GDT_Float32 : GDT_Float64, 00068 true }, 00069 // Greyscale GeoTiff scaled 0-254 00070 { "GTiff", 00071 GDT_Byte, 00072 true }, 00073 // Greyscale GeoTiff scaled 0-100 00074 { "GTiff", 00075 GDT_Byte, 00076 true }, 00077 // Greyscale BMP 00078 { "BMP", 00079 GDT_Byte, 00080 false }, 00081 // Floating HFA (Erdas Imagine Format) 00082 { "HFA", 00083 GDT_Float32, 00084 true }, 00085 // Byte HFA (Erdas Imagine Format scaled 0-100) 00086 { "HFA", 00087 GDT_Byte, 00088 true }, 00089 // Byte ASC (Arc/Info ASCII Grid Format scaled 0-100). 00090 // The GDAL code here should be AAIGrid, but this driver does not support the Create 00091 // method, so the trick is to create as HFA and then translate to AAIGrid 00092 { "HFA", 00093 GDT_Byte, 00094 true }, 00095 // Floating ASC (Arc/Info ASCII Grid floating point format). 00096 // The GDAL code here should be AAIGrid, but this driver does not support the Create 00097 // method, so the trick is to create as TIF and then translate to AAIGrid 00098 { "GTiff", 00099 (sizeof(Scalar) == 4) ? GDT_Float32 : GDT_Float64, 00100 true }, 00101 }; 00102 00103 static const GDALDataType f_type = (sizeof(Scalar) == 4) ? GDT_Float32 : GDT_Float64; 00104 00105 #ifdef MPI_FOUND 00106 int rank = 0; 00107 #endif 00108 00109 /****************************************************************/ 00110 /************************** Raster Gdal *************************/ 00111 00112 /*********************/ 00113 /*** create Raster ***/ 00114 void 00115 GdalRaster::createRaster( const string& file, int categ ) 00116 { 00117 f_file = file; 00118 f_scalefactor = 1.0; 00119 00120 // Opens read only and init the class attributes. 00121 initGdal(); 00122 open( 'r' ); 00123 00124 // set categorical status. It's not stored as part of 00125 // the gdal file's header. 00126 f_hdr.categ = categ; 00127 } 00128 00129 #ifdef MPI_FOUND 00130 void 00131 GdalRaster::createRaster( const string& output_file, const string& file, const MapFormat& format) 00132 { 00133 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 00134 00135 f_file = file; 00136 00137 // will be utilized by MPI 00138 f_output_file = output_file; 00139 00140 Scalar nv; 00141 00142 switch( format.getFormat() ) 00143 { 00144 case MapFormat::GreyBMP: 00145 f_scalefactor = 255.0; 00146 nv = 0.0; 00147 Log::instance()->debug( "GdalRaster format set to MapFormat::GreyBMP:\n" ); 00148 break; 00149 case MapFormat::GreyTiff: 00150 f_scalefactor = 254.0; 00151 nv = 255.0; 00152 Log::instance()->debug( "GdalRaster format set to MapFormat::GreyTiff:\n" ); 00153 break; 00154 case MapFormat::GreyTiff100: 00155 f_scalefactor = 100.0; 00156 nv = 127.0; //represented as 7bit so cant use 254 00157 Log::instance()->debug( "GdalRaster format set to MapFormat::GreyTiff100:\n" ); 00158 break; 00159 case MapFormat::FloatingTiff: 00160 f_scalefactor = 1.0; 00161 nv = -1.0; 00162 Log::instance()->debug( "GdalRaster format set to MapFormat::FloatingTiff:\n" ); 00163 break; 00164 case MapFormat::FloatingHFA: 00165 f_scalefactor = 1.0; 00166 nv = -1.0; 00167 Log::instance()->debug( "GdalRaster format set to MapFormat::FloatingHFA:\n" ); 00168 break; 00169 case MapFormat::ByteHFA: 00170 f_scalefactor = 100; 00171 nv = 101; 00172 Log::instance()->debug( "GdalRaster format set to MapFormat::ByteHFA:\n" ); 00173 break; 00174 case MapFormat::ByteASC: 00175 f_scalefactor = 100; 00176 nv = 101; 00177 Log::instance()->debug( "GdalRaster format set to MapFormat::ByteASC:\n" ); 00178 break; 00179 case MapFormat::FloatingASC: 00180 f_scalefactor = 1.0; 00181 nv = -9999; 00182 Log::instance()->debug( "GdalRaster format set to MapFormat::FloatingASC:\n" ); 00183 break; 00184 default: 00185 Log::instance()->error( "Unsupported output format.\n" ); 00186 throw InvalidParameterException( "Unsupported output format" ); 00187 } 00188 00189 f_hdr = Header( format.getWidth(), 00190 format.getHeight(), 00191 format.getXMin(), 00192 format.getYMin(), 00193 format.getXMax(), 00194 format.getYMax(), 00195 nv, 00196 1, 00197 0 ); 00198 00199 f_hdr.setProj( format.getProjection() ); 00200 00201 // Create a new file in disk. 00202 initGdal(); 00203 create( format.getFormat() ); 00204 } 00205 00206 #else 00207 void 00208 GdalRaster::createRaster( const string& file, const MapFormat& format) 00209 { 00210 f_file = file; 00211 00212 Scalar nv; 00213 00214 switch( format.getFormat() ) 00215 { 00216 case MapFormat::GreyBMP: 00217 f_scalefactor = 255.0; 00218 nv = 0.0; 00219 Log::instance()->debug( "GdalRaster format set to MapFormat::GreyBMP:\n" ); 00220 break; 00221 case MapFormat::GreyTiff: 00222 f_scalefactor = 254.0; 00223 nv = 255.0; 00224 Log::instance()->debug( "GdalRaster format set to MapFormat::GreyTiff:\n" ); 00225 break; 00226 case MapFormat::GreyTiff100: 00227 f_scalefactor = 100.0; 00228 nv = 127.0; //represented as 7bit so cant use 254 00229 Log::instance()->debug( "GdalRaster format set to MapFormat::GreyTiff100:\n" ); 00230 break; 00231 case MapFormat::FloatingTiff: 00232 f_scalefactor = 1.0; 00233 nv = -1.0; 00234 Log::instance()->debug( "GdalRaster format set to MapFormat::FloatingTiff:\n" ); 00235 break; 00236 case MapFormat::FloatingHFA: 00237 f_scalefactor = 1.0; 00238 nv = -1.0; 00239 Log::instance()->debug( "GdalRaster format set to MapFormat::FloatingHFA:\n" ); 00240 break; 00241 case MapFormat::ByteHFA: 00242 f_scalefactor = 100; 00243 nv = 101; 00244 Log::instance()->debug( "GdalRaster format set to MapFormat::ByteHFA:\n" ); 00245 break; 00246 case MapFormat::ByteASC: 00247 f_scalefactor = 100; 00248 nv = 101; 00249 Log::instance()->debug( "GdalRaster format set to MapFormat::ByteASC:\n" ); 00250 break; 00251 case MapFormat::FloatingASC: 00252 f_scalefactor = 1.0; 00253 nv = -9999; 00254 Log::instance()->debug( "GdalRaster format set to MapFormat::FloatingASC:\n" ); 00255 break; 00256 default: 00257 Log::instance()->error( "Unsupported output format.\n" ); 00258 throw InvalidParameterException( "Unsupported output format" ); 00259 } 00260 00261 f_hdr = Header( format.getWidth(), 00262 format.getHeight(), 00263 format.getXMin(), 00264 format.getYMin(), 00265 format.getXMax(), 00266 format.getYMax(), 00267 nv, 00268 1, 00269 0 ); 00270 00271 f_hdr.setProj( format.getProjection() ); 00272 00273 // Create a new file in disk. 00274 initGdal(); 00275 create( format.getFormat() ); 00276 } 00277 #endif 00278 00279 /*****************/ 00280 /*** destructor ***/ 00281 GdalRaster::~GdalRaster() 00282 { 00283 #ifdef MPI_FOUND 00284 if (((rank != 0) && (strcmp (f_file.c_str(), f_output_file.c_str()) != 0)) || (rank == 0)){ 00285 #endif 00286 // Save the last line read, if needed. 00287 saveRow(); 00288 00289 if ( f_data ) { 00290 00291 delete[] f_data; 00292 } 00293 00294 GDALClose( f_ds ); 00295 #ifdef MPI_FOUND 00296 } 00297 #endif 00298 } 00299 00300 00301 /************/ 00302 /*** read ***/ 00303 void 00304 GdalRaster::read( Scalar *buf, int frow, int nrow ) 00305 { 00306 // Header's data 00307 int nband = f_hdr.nband; 00308 int size = f_size * nrow; 00309 00310 // Read each band 00311 for ( int i = 1; i <= nband; i++, buf += size ) { 00312 00313 // Gets the i-th band. 00314 GDALRasterBand *band = f_ds->GetRasterBand( i ); 00315 00316 // Fill 'buf' with new data 00317 int ret = band->RasterIO( GF_Read, 0, frow, f_size, nrow, buf, f_size, nrow, f_type, 0, 0 ); 00318 00319 if ( ret == CE_Failure ) { 00320 00321 // Fill 'buf' with nodata 00322 for ( int j = 0 ; j < f_size; j++ ) { 00323 00324 buf[j] = f_hdr.noval; 00325 } 00326 } 00327 } 00328 } 00329 00330 00331 /*************/ 00332 /*** write ***/ 00333 void 00334 GdalRaster::write( Scalar *buf, int frow, int nrow ) 00335 { 00336 // Header's data. 00337 int nband = f_hdr.nband; 00338 int size = f_size * nrow; 00339 00340 // Write each band in the file. 00341 for ( int i = 1; i <= nband; i++, buf += size ) { 00342 00343 // Gets the i-th band. 00344 GDALRasterBand *band = f_ds->GetRasterBand( i ); 00345 00346 // Write the buffer's data in the file. 00347 int ret = band->RasterIO( GF_Write, 0, frow, f_size, nrow, buf, f_size, nrow, f_type, 0, 0 ); 00348 00349 if ( ret == CE_Failure ) 00350 { 00351 Log::instance()->error( "Unable to write to file %s\n", f_file.c_str()); 00352 throw FileIOException( "Unable to write to file " + f_file, f_file ); 00353 } 00354 } 00355 } 00356 00357 /************/ 00358 /*** open ***/ 00359 void 00360 GdalRaster::open( char mode ) 00361 { 00362 // Mode: write or read. 00363 GDALAccess gmod = (mode == 'w') ? GA_Update : GA_ReadOnly; 00364 00365 // Opens the file. 00366 f_ds = (GDALDataset *)GDALOpenShared( f_file.c_str(), gmod ); 00367 00368 if ( ! f_ds ) { 00369 00370 Log::instance()->error( "Unable to open file %s\n", f_file.c_str()); 00371 throw FileIOException( "Unable to open file " + f_file, f_file ); 00372 } 00373 00374 // Number of bands (channels). 00375 f_hdr.nband = f_ds->GetRasterCount(); 00376 00377 // Projection. 00378 const char * projwkt = (char *)f_ds->GetProjectionRef(); 00379 00380 f_hdr.setProj( projwkt ); 00381 00382 if ( ! f_hdr.hasProj() ) { 00383 00384 Log::instance()->warn( "The raster %s is not georeferenced. Assuming LatLong WGS84\n", f_file.c_str() ); 00385 f_hdr.setProj( GeoTransform::getDefaultCS() ); 00386 } 00387 00388 // Assumes that all bands have the same georeference 00389 // characteristics, ie, the same header. 00390 00391 GDALRasterBand *band = f_ds->GetRasterBand( 1 ); 00392 00393 int xd = f_hdr.xdim = band->GetXSize(); 00394 int yd = f_hdr.ydim = band->GetYSize(); 00395 00396 f_hdr.noval = band->GetNoDataValue(); 00397 00398 // Map limits. 00399 double cf[6]; 00400 f_ds->GetGeoTransform( cf ); 00401 f_hdr.xmin = (Scalar) cf[0]; 00402 f_hdr.xmax = (Scalar) (cf[0] + xd * cf[1] + yd * cf[2]); 00403 f_hdr.ymax = (Scalar) cf[3]; 00404 f_hdr.ymin = (Scalar) (cf[3] + xd * cf[4] + yd * cf[5]); 00405 00406 //Log::instance()->debug( "Raster boundaries: xmin=%f, xmax=%f, ymin=%f, ymax=%f \n", f_hdr.xmin, f_hdr.xmax, f_hdr.ymin, f_hdr.ymax ); 00407 00408 f_hdr.calculateCell(); 00409 00410 for ( int i=0; i<6; ++i ) { 00411 00412 f_hdr.gt[i] = cf[i]; 00413 } 00414 00415 // Minimum and maximum values. 00416 f_hdr.vmin = band->GetMinimum( &f_hdr.minmax ); 00417 00418 if ( f_hdr.minmax ) { 00419 00420 f_hdr.vmax = band->GetMaximum(); 00421 } 00422 00423 // Assumes grid (in place of 'pixel'). 00424 // todo: how can I read this with GDAL? 00425 f_hdr.grid = 0; 00426 00427 // Create warped data set when coordinate system is different 00428 if ( ! GeoTransform::compareCoordSystemStrings( projwkt, GeoTransform::getDefaultCS() ) ) { 00429 00430 GDALWarpOptions *warpOptions = GDALCreateWarpOptions(); 00431 00432 f_warped_ds = (GDALDataset *)GDALAutoCreateWarpedVRT( f_ds, projwkt, GeoTransform::getDefaultCS(), GRA_NearestNeighbour, 0.0, warpOptions ); 00433 } 00434 00435 // Initialize the Buffer 00436 initBuffer(); 00437 } 00438 00439 /**************/ 00440 /*** create ***/ 00441 void 00442 GdalRaster::create( int format ) 00443 { 00444 // Store format for future reference (used in method "finish") 00445 f_format = format; 00446 00447 GDALDriver *poDriver; 00448 char **papszMetadata; 00449 00450 char const *fmt = Formats[ format ].GDalDriverName; 00451 00452 // Added by Tim to see if the chosen format supports the gdal create method 00453 poDriver = GetGDALDriverManager()->GetDriverByName(fmt); 00454 00455 if ( poDriver == NULL ) { 00456 00457 std::string msg = "Unable to load GDAL driver " + string(fmt) + " for file " + f_file; 00458 00459 Log::instance()->error( msg.c_str() ); 00460 00461 throw FileIOException( "Unable to load GDAL driver " + string(fmt), f_file ); 00462 } 00463 00464 papszMetadata = poDriver->GetMetadata(); 00465 00466 if ( ! CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) ) { 00467 00468 Log::instance()->error( "Driver %s, format %s does not support Create() method.\n", 00469 poDriver->GetDescription(), 00470 fmt ); 00471 00472 throw FileIOException( "Driver " + string(fmt) + " does not support Create() method", f_file ); 00473 } 00474 00475 // Read the parameters in 'hdr' used to create the file. 00476 // The other parameters come from the copied header. 00477 f_hdr.nband = 1; 00478 00479 // Create the file. 00480 if ( format == MapFormat::FloatingHFA || format == MapFormat::ByteHFA ) { 00481 00482 // Note: HFA (erdas imagine) format does not support nodata before GDAL 1.5 00483 // see http://www.gdal.org/gdal_tutorial.html for options examples 00484 char **papszOptions = NULL; 00485 00486 papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "YES" ); 00487 #ifdef MPI_FOUND 00488 if (rank ==0){ 00489 #endif 00490 f_ds = poDriver->Create( f_file.c_str(), 00491 f_hdr.xdim, f_hdr.ydim, 00492 f_hdr.nband, 00493 Formats[format].dataType, 00494 papszOptions ); 00495 #ifdef MPI_FOUND 00496 } 00497 #endif 00498 CSLDestroy( papszOptions ); 00499 } 00500 //ArcMap needs a LZW compression license to read files compressed 00501 //like this so you may need to comment out this code if you want 00502 //to open the generated maps without having the license. Compression 00503 //is enabled for now because it offers significant size reduction. 00504 else if (format == MapFormat::GreyTiff100) 00505 { 00506 //lzw compression and represent each pixel with 7bits only 00507 char **papszOptions = NULL; 00508 papszOptions = CSLSetNameValue( papszOptions, "NBITS", "7" ); 00509 papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "LZW" ); 00510 #ifdef MPI_FOUND 00511 if (rank ==0){ 00512 #endif 00513 f_ds = poDriver->Create( f_file.c_str(), 00514 f_hdr.xdim, f_hdr.ydim, 00515 f_hdr.nband, 00516 Formats[format].dataType, //data type 00517 papszOptions ); //opt parameters 00518 #ifdef MPI_FOUND 00519 } 00520 #endif 00521 CSLDestroy( papszOptions ); 00522 } 00523 //Create temporary ByteHFA file with a different name (original file name + .tmp) 00524 //It will be converted to ASC in the finish method 00525 else if (format==MapFormat::ByteASC) 00526 { 00527 std::string temp_file = f_file + ".tmp"; 00528 00529 #ifdef MPI_FOUND 00530 if (rank ==0){ 00531 #endif 00532 f_ds = poDriver->Create( temp_file.c_str(), 00533 f_hdr.xdim, f_hdr.ydim, 00534 f_hdr.nband, 00535 Formats[format].dataType, 00536 NULL ); 00537 #ifdef MPI_FOUND 00538 } 00539 #endif 00540 } 00541 //Create temporary GeoTiff file with a different name (original file name + .tmp) 00542 //It will be converted to ASC in the finish method 00543 else if (format==MapFormat::FloatingASC) 00544 { 00545 std::string temp_file = f_file + ".tmp"; 00546 00547 #ifdef MPI_FOUND 00548 if (rank ==0){ 00549 #endif 00550 f_ds = poDriver->Create( temp_file.c_str(), 00551 f_hdr.xdim, f_hdr.ydim, 00552 f_hdr.nband, 00553 Formats[format].dataType, 00554 NULL ); 00555 #ifdef MPI_FOUND 00556 } 00557 #endif 00558 } 00559 else { 00560 #ifdef MPI_FOUND 00561 if (rank ==0){ 00562 #endif 00563 //uncompressed 00564 f_ds = poDriver->Create( f_file.c_str(), 00565 f_hdr.xdim, f_hdr.ydim, 00566 f_hdr.nband, 00567 Formats[format].dataType, 00568 NULL ); 00569 #ifdef MPI_FOUND 00570 } 00571 #endif 00572 } 00573 00574 #ifdef MPI_FOUND 00575 if (rank ==0) { 00576 #endif 00577 if ( ! f_ds ) { 00578 00579 Log::instance()->warn( "Unable to create file %s.\n",f_file.c_str() ); 00580 throw FileIOException( "Unable to create file " + f_file, f_file ); 00581 } 00582 00583 if ( Formats[format].hasMeta ) { 00584 00585 // Metadata 00586 00587 // Limits (without rotations). 00588 f_ds->SetGeoTransform( f_hdr.gt ); 00589 00590 // Projection. 00591 f_ds->SetProjection( f_hdr.proj.c_str() ); 00592 00593 // Sets the "nodata" value in all bands. 00594 int ret = f_ds->GetRasterBand(1)->SetNoDataValue( f_hdr.noval ); 00595 00596 if ( ret == CE_Failure ) { 00597 00598 Log::instance()->warn( "Raster %s (%s) does not support nodata value assignment. Nodata values will correspond to %f anyway, but this will not be stored as part of the raster metadata.\n", f_file.c_str(), fmt, f_hdr.noval ); 00599 } 00600 } 00601 #ifdef MPI_FOUND 00602 } 00603 #endif 00604 // Initialize the Buffer 00605 initBuffer(); 00606 } 00607 00608 /*******************/ 00609 /*** init Buffer ***/ 00610 void 00611 GdalRaster::initBuffer() 00612 { 00613 // Initialize the buffer. 00614 f_size = f_hdr.xdim; 00615 00616 if ( f_data ) { 00617 00618 delete[] f_data; 00619 } 00620 00621 f_data = new Scalar[ f_size * f_hdr.nband ]; 00622 00623 f_changed = 0; 00624 f_currentRow = -1; 00625 } 00626 00627 /*****************/ 00628 /*** init Gdal ***/ 00629 void 00630 GdalRaster::initGdal() 00631 { 00632 static int first = 1; 00633 00634 if ( first ) { 00635 00636 first = 0; 00637 GDALAllRegister(); 00638 } 00639 } 00640 00641 /************/ 00642 /*** iput ***/ 00643 int 00644 GdalRaster::iput( int x, int y, Scalar val ) 00645 { 00646 // Be sure that 'y' line is in the buffer 'f_data'. 00647 loadRow( y, true ); 00648 00649 // Put values in the first band of (x,y) position. 00650 f_data[x] = val; 00651 00652 // Indicates the line 'f_currentRow' has changed. 00653 f_changed = 1; 00654 00655 return 1; 00656 } 00657 00658 /************/ 00659 /*** iget ***/ 00660 int 00661 GdalRaster::iget( int x, int y, Scalar *val ) 00662 { 00663 // Be sure that 'y' line is in the buffer 'f_data'. 00664 loadRow( y ); 00665 00666 // Get all band's values. 00667 Scalar *pv = val; 00668 int nband = f_hdr.nband; 00669 for ( int i = 0; i < nband; i++, x += f_size ) { 00670 00671 if ( (*pv++ = f_data[x]) == f_hdr.noval ) { 00672 00673 return 0; 00674 } 00675 } 00676 00677 return 1; 00678 } 00679 00680 /***********/ 00681 /*** get ***/ 00682 int 00683 GdalRaster::get( Coord px, Coord py, Scalar *val ) 00684 { 00685 Scalar *pv = val; 00686 for ( int i = 0; i < f_hdr.nband; i++ ) { 00687 00688 *pv++ = f_hdr.noval; 00689 } 00690 00691 pair<int,int> xy = f_hdr.convertLonLat2XY( px, py ); 00692 int x = xy.first; 00693 int y = xy.second; 00694 00695 // If the point is out of range, returns 0. 00696 if ( x < 0 || x >= f_hdr.xdim || y < 0 || y >= f_hdr.ydim ) { 00697 00698 //Log::instance()->debug( "Raster::get() Pixel (%d,%d) is not in extent\n",x,y); 00699 return 0; 00700 } 00701 00702 // 'iget()' detects if the coordinate has or not valid 00703 // information (noval); 00704 return iget( x, y, val ); 00705 } 00706 00707 /***********/ 00708 /*** put ***/ 00709 int 00710 GdalRaster::put( Coord px, Coord py, Scalar val ) 00711 { 00712 pair<int,int> xy = f_hdr.convertLonLat2XY( px, py ); 00713 int x = xy.first; 00714 int y = xy.second; 00715 00716 if ( x < 0 || x >= f_hdr.xdim || y < 0 || y >= f_hdr.ydim ) { 00717 00718 Log::instance()->warn( "Coordinate (%f, %f) corresponds to a cell (%d, %d) outside the raster boundaries. It will be ignored.\n", px, py, x, y ); 00719 return 0; 00720 } 00721 00722 return iput( x, y, f_scalefactor*val ); 00723 } 00724 00725 /***********/ 00726 /*** put ***/ 00727 int 00728 GdalRaster::put( Coord px, Coord py ) 00729 { 00730 pair<int,int> xy = f_hdr.convertLonLat2XY( px, py ); 00731 int x = xy.first; 00732 int y = xy.second; 00733 00734 Scalar val = f_hdr.noval; 00735 00736 if ( x < 0 || x >= f_hdr.xdim || y < 0 || y >= f_hdr.ydim ) { 00737 00738 Log::instance()->warn( "Coordinate (%f, %f) corresponds to a cell (%d, %d) outside the raster boundaries. It will be ignored.\n", px, py, x, y ); 00739 return 0; 00740 } 00741 00742 return iput( x, y, val ); 00743 } 00744 00745 /*******************/ 00746 /*** get Min Max ***/ 00747 int 00748 GdalRaster::getMinMax( Scalar *min, Scalar *max ) 00749 { 00750 if ( ! calcMinMax() ) { 00751 00752 return 0; 00753 } 00754 00755 *min = f_hdr.vmin; 00756 00757 *max = f_hdr.vmax; 00758 00759 return 1; 00760 } 00761 00762 /********************/ 00763 /*** calc Min Max ***/ 00764 int 00765 GdalRaster::calcMinMax( int band ) 00766 { 00767 if ( f_hdr.minmax ) { 00768 00769 return 1; 00770 } 00771 00772 Scalar *bands; 00773 bands = new Scalar[ f_hdr.nband ]; 00774 Scalar *val = bands + band; 00775 Scalar min=0; 00776 Scalar max=0; 00777 00778 bool initialized = false; 00779 00780 // Scan all map values. 00781 for ( int y = 0; y < f_hdr.ydim; y++ ) { 00782 00783 for ( int x = 0; x < f_hdr.xdim; x++ ) { 00784 00785 if ( iget( x, y, bands ) ) { 00786 00787 if ( !initialized ) { 00788 00789 initialized = true; 00790 min = max = *val; 00791 } 00792 00793 if ( min > *val ) { 00794 00795 min = *val; 00796 } 00797 else if ( max < *val ) { 00798 00799 max = *val; 00800 } 00801 } 00802 } 00803 } 00804 00805 delete[] bands; 00806 00807 if ( ! initialized ) { 00808 00809 return 0; 00810 } 00811 00812 // This is only logically const. Here we cast away const to cache the min/max 00813 Header *f_hdrNoConst = const_cast<Header*>(&f_hdr); 00814 00815 f_hdrNoConst->minmax = 1; 00816 f_hdrNoConst->vmin = min; 00817 f_hdrNoConst->vmax = max; 00818 00819 return 1; 00820 } 00821 00822 /****************/ 00823 /*** load Row ***/ 00824 void 00825 GdalRaster::loadRow( int row, bool writeOperation ) 00826 { 00827 // If the line is already read, return. 00828 if ( row == f_currentRow ) { 00829 00830 return; 00831 } 00832 00833 saveRow(); 00834 00835 // Update f_data 00836 if ( writeOperation ) { 00837 00838 // Just reset f_data with nodata 00839 for ( int i = 0 ; i < f_size; i++ ) { 00840 00841 f_data[i] = f_hdr.noval; 00842 } 00843 } 00844 else { 00845 00846 // Read from file 00847 read( f_data, row, 1 ); 00848 } 00849 00850 f_currentRow = row; 00851 f_changed = 0; 00852 } 00853 00854 00855 /****************/ 00856 /*** save Row ***/ 00857 void 00858 GdalRaster::saveRow() 00859 { 00860 if ( ! f_changed || f_currentRow < 0 ) { 00861 00862 return; 00863 } 00864 00865 write( f_data, f_currentRow, 1 ); 00866 00867 f_changed = 0; 00868 } 00869 00870 /**************/ 00871 /*** finish ***/ 00872 void 00873 GdalRaster::finish() 00874 { 00875 // Save the last line read, if needed. 00876 saveRow(); 00877 00878 if ( f_format == MapFormat::ByteASC || f_format == MapFormat::FloatingASC ) 00879 { 00880 // Temporary file name 00881 std::string temp_file = f_file + ".tmp"; 00882 00883 // Get ArcInfo/ASC Grid driver 00884 GDALDriver *hDriver; 00885 hDriver = GetGDALDriverManager()->GetDriverByName("AAIGrid"); 00886 00887 if ( hDriver == NULL ) { 00888 00889 std::string msg = "Unable to load AAIGrid GDAL driver"; 00890 00891 Log::instance()->error( msg.c_str() ); 00892 throw FileIOException( msg.c_str(), f_file ); 00893 } 00894 00895 // Convert temporary raster to ArcInfo/ASC Grid 00896 GDALDataset * new_ds = hDriver->CreateCopy( f_file.c_str(), f_ds, FALSE, NULL, NULL, NULL ); 00897 00898 if ( ! new_ds ) { 00899 00900 Log::instance()->warn( "Unable to create raster copy %s.\n",f_file.c_str() ); 00901 throw FileIOException( "Unable to create raster copy " + f_file, f_file ); 00902 } 00903 00904 // Delete temporary ByteHFA raster 00905 #if GDAL_VERSION_MAJOR > 1 || ( GDAL_VERSION_MAJOR == 1 && GDAL_VERSION_MINOR >= 5 ) 00906 00907 delete f_ds; 00908 00909 if ( GDALDriver::QuietDelete( temp_file.c_str() ) == CE_Failure ) 00910 { 00911 Log::instance()->warn( "Could not delete temporary file %s", temp_file.c_str() ); 00912 } 00913 #else 00914 GDALDriver * temp_driver = f_ds->GetDriver(); 00915 00916 if ( temp_driver->Delete( temp_file.c_str() ) == CE_Failure ) 00917 { 00918 Log::instance()->warn( "Could not delete temporary file %s", temp_file.c_str() ); 00919 } 00920 00921 delete f_ds; 00922 #endif 00923 00924 f_ds = new_ds; 00925 } 00926 } 00927 00928 /*********************/ 00929 /*** delete Raster ***/ 00930 int 00931 GdalRaster::deleteRaster() 00932 { 00933 GDALDriver * driver = f_ds->GetDriver(); 00934 00935 int ret = driver->Delete( f_file.c_str() ); 00936 00937 if ( ret == CE_Failure ) { 00938 00939 return 0; 00940 } 00941 00942 f_ds = 0; 00943 00944 return 1; 00945 } 00946 00947 /*****************************/ 00948 /*** getExtentInStandardCs ***/ 00949 int 00950 GdalRaster::getExtentInStandardCs( Coord *xmin, Coord *ymin, Coord *xmax, Coord *ymax ) 00951 { 00952 if ( f_warped_ds == 0 ) { 00953 00954 return 0; 00955 } 00956 00957 GDALRasterBand *band = f_warped_ds->GetRasterBand( 1 ); 00958 00959 int xd = band->GetXSize(); 00960 int yd = band->GetYSize(); 00961 00962 double cf[6]; 00963 f_warped_ds->GetGeoTransform( cf ); 00964 00965 *xmin = cf[0]; 00966 *ymin = cf[3] + xd * cf[4] + yd * cf[5]; 00967 *xmax = cf[0] + xd * cf[1] + yd * cf[2]; 00968 *ymax = cf[3]; 00969 00970 return 1; 00971 }