openModeller
Version 1.4.0
|
00001 00030 #include <openmodeller/env_io/TerralibRaster.hh> 00031 00032 #include <openmodeller/TeDatabaseManager.hh> 00033 #include <openmodeller/TeStringParser.hh> 00034 #include <TeRaster.h> 00035 #include <TeDatabase.h> 00036 #include <TeLayer.h> 00037 00038 #include <openmodeller/env_io/GeoTransform.hh> 00039 #include <openmodeller/Log.hh> 00040 #include <openmodeller/Exceptions.hh> 00041 #include <openmodeller/MapFormat.hh> 00042 00043 #include <TeProjection.h> 00044 00045 #include <string> 00046 using std::string; 00047 00048 #include <vector> 00049 using std::vector; 00050 00051 #include <utility> 00052 using std::pair; 00053 00054 #include <stdio.h> 00055 /****************************************************************/ 00056 /************************** Te Raster ***************************/ 00060 Raster* 00061 TerralibRaster::CreateRasterCallback() 00062 { 00063 return new TerralibRaster(); 00064 } 00065 00071 // Open an existing file -- read only. 00072 void 00073 TerralibRaster::createRaster( const string& str, int categ ) 00074 { 00075 te_str_parser_ = new TeStringParser(); 00076 te_str_parser_->str_ = str; 00077 f_scalefactor = 1.0; 00078 f_hdr.minmax = 0; 00079 00080 // Open an existing Raster 00081 openTeRaster(); 00082 00083 if ( raster_ ) 00084 { 00085 00086 raster_->init(); 00087 params_ = &raster_->params(); 00088 00089 if (params_->status_ != TeRasterParams::TeReadyToRead) 00090 { 00091 std::string msg = "TerralibRaster::createRaster - Raster cannot be opened: "; 00092 msg += raster_->errorMessage().c_str(); 00093 Log::instance()->error( msg.c_str() ); 00094 throw RasterException( msg ); 00095 } 00096 else 00097 { 00098 // Number of bands (channels). 00099 f_hdr.nband = raster_->nBands(); 00100 00101 // Projection. 00102 //string teste = TeGetWKTFromTeProjection( raster_->projection() ); 00103 f_hdr.setProj( TeGetWKTFromTeProjection( raster_->projection() ) ); 00104 if ( ! f_hdr.hasProj() ) 00105 { 00106 Log::instance()->warn( "Raster %s is not georeferenced. Assuming LatLong WGS84\n", f_file.c_str() ); 00107 f_hdr.setProj( GeoTransform::getDefaultCS() ); 00108 } 00109 00110 // Assumes that all bands have the same georeference 00111 // characteristics, ie, the same header. 00112 f_hdr.xdim = params_->ncols_; 00113 f_hdr.ydim = params_->nlines_; 00114 00115 f_hdr.noval = params_->dummy_[0]; 00116 00117 // Map limits. 00118 f_hdr.xmin = (Scalar) params_->boundingBox().x1_; 00119 f_hdr.xmax = (Scalar) params_->boundingBox().x2_; 00120 f_hdr.ymax = (Scalar) params_->boundingBox().y2_; 00121 f_hdr.ymin = (Scalar) params_->boundingBox().y1_; 00122 00123 f_hdr.calculateCell(); 00124 00125 f_hdr.gt[0] = params_->boundingBox().x1_; // top left x 00126 f_hdr.gt[1] = params_->resx_; // w-e pixel resolution 00127 f_hdr.gt[2] = 0; // rotation, 0 if image is "north up" 00128 f_hdr.gt[3] = params_->boundingBox().y2_; // top left y 00129 f_hdr.gt[4] = 0; // rotation, 0 if image is "north up" 00130 f_hdr.gt[5] = - params_->resy_; // n-s pixel resolution 00131 00132 // Minimum and maximum values. 00133 f_hdr.min = params_->vmin_[0]; 00134 f_hdr.max = params_->vmax_[0]; 00135 00136 // If not zero 'min' and 'max' are valid values. 00137 f_hdr.minmax = 1; //true 00138 00139 // Assumes grid (in place of 'pixel'). 00140 // todo: how can I read this with TerraLib? 00141 f_hdr.grid = 0; 00142 00143 // set categorical status. It's not stored as part of 00144 // the gdal file's header. 00145 f_hdr.categ = categ; 00146 } 00147 } 00148 } 00149 00155 void 00156 TerralibRaster::createRaster( const std::string& str, const MapFormat& format ) 00157 { 00158 te_str_parser_ = new TeStringParser(); 00159 te_str_parser_->str_ = str; 00160 00161 Scalar nv; 00162 00163 TeDataType rtype; 00164 switch( format.getFormat() ) 00165 { 00166 case MapFormat::GreyBMP: 00167 f_scalefactor = 255.0; 00168 nv = 0.0; 00169 rtype = TeUNSIGNEDCHAR; 00170 break; 00171 case MapFormat::GreyTiff: 00172 f_scalefactor = 254.0; 00173 nv = 255.0; 00174 rtype = TeUNSIGNEDCHAR; 00175 break; 00176 case MapFormat::FloatingTiff: 00177 f_scalefactor = 1.0; 00178 nv = -1.0; 00179 rtype = TeFLOAT; 00180 break; 00181 default: 00182 Log::instance()->error( "Unsupported output format.\n" ); 00183 throw InvalidParameterException( "Unsupported output format" ); 00184 } 00185 00186 f_hdr = Header( format.getWidth(), 00187 format.getHeight(), 00188 format.getXMin(), 00189 format.getYMin(), 00190 format.getXMax(), 00191 format.getYMax(), 00192 nv, 00193 1, 00194 0 ); 00195 00196 f_hdr.setProj( format.getProjection() ); 00197 00198 // Create a new TeRaster. 00199 createTeRaster(); 00200 00201 if ( raster_ ) 00202 { 00203 params_ = &raster_->params(); 00204 00205 params_->nBands( 1 ); 00206 params_->setDataType( rtype ); 00207 params_->setDummy(nv, 0); 00208 00209 params_->boundingBoxLinesColumns(format.getXMin(), format.getYMin(), 00210 format.getXMax(), format.getYMax(), 00211 format.getHeight(), format.getWidth()); 00212 00213 // Block Size. 00214 /*if(format.getHeight() > 128 || format.getWidth() > 128) 00215 { 00216 params_->blockHeight_ = 128; 00217 params_->blockWidth_ = 128; 00218 } 00219 else 00220 { 00221 params_->blockHeight_ = format.getHeight(); 00222 params_->blockWidth_ = format.getWidth(); 00223 }*/ 00224 params_->blockHeight_ = 1; 00225 params_->blockWidth_ = format.getWidth(); 00226 00227 raster_->init(); 00228 00229 if (params_->status_ != TeRasterParams::TeReadyToWrite) 00230 { 00231 Log::instance()->warn( "TerraLib raster not ready: %s.\n", raster_->errorMessage() ); 00232 } 00233 } 00234 } 00235 00239 TerralibRaster::~TerralibRaster() 00240 { 00241 if( layer_ ) 00242 { 00243 layer_->addRasterGeometry( raster_ ); 00244 } 00245 00246 delete te_str_parser_; 00247 } 00248 00253 int 00254 TerralibRaster::get( Coord px, Coord py, Scalar *val ) 00255 { 00256 pair<int,int> xy = f_hdr.convertLonLat2XY( px, py ); 00257 int x = xy.first; 00258 int y = xy.second; 00259 00260 // If the point is out of range, returns 0. 00261 if ( x < 0 || x >= f_hdr.xdim || y < 0 || y >= f_hdr.ydim ) 00262 { 00263 //Log::instance()->debug( "Raster::get() Pixel (%d,%d) is not in extent\n",x,y); 00264 return 0; 00265 } 00266 00267 if( raster_->getElement(x, y, *val , 0) ) 00268 return 1; 00269 00270 return 0; 00271 }; 00272 00278 int 00279 TerralibRaster::put( Coord px, Coord py, Scalar val ) 00280 { 00281 pair<int,int> xy = f_hdr.convertLonLat2XY( px, py ); 00282 int x = xy.first; 00283 int y = xy.second; 00284 00285 if ( x < 0 || x >= f_hdr.xdim || y < 0 || y >= f_hdr.ydim ) 00286 return 0; 00287 00288 if( raster_->setElement(x, y, f_scalefactor*val, 0) ) 00289 return 1; 00290 00291 return 0; 00292 } 00293 00299 int 00300 TerralibRaster::put( Coord px, Coord py ) 00301 { 00302 pair<int,int> xy = f_hdr.convertLonLat2XY( px, py ); 00303 int x = xy.first; 00304 int y = xy.second; 00305 00306 if ( x < 0 || x >= f_hdr.xdim || y < 0 || y >= f_hdr.ydim ) 00307 return 0; 00308 00309 if( raster_->setElement(x, y, f_hdr.noval, 0) ) 00310 return 1; 00311 00312 return 0; 00313 } 00314 00318 void 00319 TerralibRaster::openTeRaster() 00320 { 00321 if( te_str_parser_->parse() ) 00322 { 00323 db_ = TeDatabaseManager::instance().create( *te_str_parser_ ); 00324 00325 if ( !db_->isConnected() ) 00326 { 00327 //delete db_; 00328 std::string msg = "TerralibRaster::openTeRaster - Cannot connect to database: "; 00329 msg += db_->errorMessage().c_str(); 00330 Log::instance()->error( msg.c_str() ); 00331 throw RasterException( msg ); 00332 } 00333 else 00334 { 00335 if (db_->layerExist( te_str_parser_->layerName_ )) 00336 { 00337 layer_ = new TeLayer(te_str_parser_->layerName_, db_); 00338 raster_ = layer_->raster(); 00339 } 00340 } 00341 } 00342 else 00343 { 00344 // Disk file. 00345 raster_ = new TeRaster( te_str_parser_->str_, 'r'); 00346 } 00347 } 00348 00352 void 00353 TerralibRaster::createTeRaster() 00354 { 00355 if( te_str_parser_->parse() ) 00356 { 00357 db_ = TeDatabaseManager::instance().create( *te_str_parser_ ); 00358 00359 if ( !db_->isConnected() ) 00360 { 00361 //delete db_; 00362 std::string msg = "TerralibRaster::createTeRaster - Cannot connect to database: "; 00363 msg += db_->errorMessage().c_str(); 00364 Log::instance()->error( msg.c_str() ); 00365 throw RasterException( msg ); 00366 } 00367 else 00368 { 00369 TeProjection* proj = TeGetTeProjectionFromWKT( f_hdr.proj ); 00370 00371 layer_ = new TeLayer(te_str_parser_->layerName_, db_, proj); 00372 00373 // create a raster geometry in a TerraLib database 00374 TeRasterParams params; 00375 params.fileName_ = te_str_parser_->layerName_; 00376 params.mode_ = 'c'; 00377 // parameters specific of the database decoder 00378 params.decoderIdentifier_ = "DB"; // a database decoder 00379 params.database_ = db_; // pointer to the database 00380 // we do not intent to mosaic any raster data to this 00381 params.tiling_type_ = TeRasterParams::TeNoExpansible; 00382 00383 // the photometric interpretation of the raster 00384 params.setPhotometric(TeRasterParams::TeMultiBand); 00385 00386 raster_ = new TeRaster( params ); 00387 } 00388 } 00389 else 00390 { 00391 // Disk file. 00392 raster_ = new TeRaster( te_str_parser_->str_, 'c'); 00393 } 00394 } 00395 00396 /*******************/ 00397 /*** get Min Max ***/ 00398 int 00399 TerralibRaster::getMinMax( Scalar *min, Scalar *max ) 00400 { 00401 *min = f_hdr.min; 00402 00403 *max = f_hdr.max; 00404 00405 return 1; 00406 } 00407 00408 00409 /*********************/ 00410 /*** delete Raster ***/ 00411 int 00412 TerralibRaster::deleteRaster() 00413 { 00414 db_ = TeDatabaseManager::instance().create( *te_str_parser_ ); 00415 00416 if ( !db_->isConnected() ) 00417 { 00418 //delete db_; 00419 std::string msg = "TerralibRaster::openTeRaster - Cannot connect to database: "; 00420 msg += db_->errorMessage().c_str(); 00421 Log::instance()->error( msg.c_str() ); 00422 throw RasterException( msg ); 00423 return 0; 00424 } 00425 else 00426 { 00427 if ( db_->layerExist( layer_->name() ) ) 00428 { 00429 db_->deleteLayer( layer_->id() ); 00430 } 00431 00432 if( raster_ ) 00433 { 00434 delete raster_; 00435 } 00436 00437 if( params_ ) 00438 { 00439 delete params_; 00440 } 00441 00442 if( te_str_parser_ ) 00443 { 00444 delete te_str_parser_; 00445 } 00446 00447 return 1; 00448 } 00449 00450 // If raster is a file disk. 00451 if( raster_ ) 00452 { 00453 delete raster_; 00454 } 00455 00456 if( params_ ) 00457 { 00458 delete params_; 00459 } 00460 00461 if( te_str_parser_ ) 00462 { 00463 delete te_str_parser_; 00464 } 00465 00466 return 1; 00467 }