openModeller
Version 1.4.0
|
00001 00028 #include "aquamaps.hh" 00029 00030 #include <openmodeller/Configuration.hh> 00031 #include <openmodeller/Exceptions.hh> 00032 #include <openmodeller/os_specific.hh> 00033 00034 #include <string.h> 00035 #include <stdio.h> 00036 #include <math.h> 00037 #include <sstream> 00038 #include <fstream> 00039 #include <sqlite3.h> 00040 00041 // Algorithm is included for std::min and std::max. 00042 #include <algorithm> 00043 00044 using namespace std; 00045 00046 00047 /****************************************************************/ 00048 /********************** Algorithm's Metadata ********************/ 00049 00050 #define NUM_PARAM 7 00051 00052 /*************************************/ 00053 /*** Algorithm parameters metadata ***/ 00054 00055 #define PARAM_USE_SURFACE_LAYERS "UseSurfaceLayers" 00056 #define PARAM_USE_DEPTH_RANGE "UseDepthRange" 00057 #define PARAM_USE_ICE_CONCENTRATION "UseIceConcentration" 00058 #define PARAM_USE_DISTANCE_TO_LAND "UseDistanceToLand" 00059 #define PARAM_USE_PRIMARY_PRODUCTION "UsePrimaryProduction" 00060 #define PARAM_USE_SALINITY "UseSalinity" 00061 #define PARAM_USE_TEMPERATURE "UseTemperature" 00062 00063 static AlgParamMetadata parameters[NUM_PARAM] = { // Parameters 00064 { 00065 PARAM_USE_SURFACE_LAYERS, // Id 00066 "Use surface layers (only for temperature and salinity)", // Name 00067 Integer, // Type 00068 "Use surface layers (1=yes, 0=no, -1=let the algorithm decide)", // Overview 00069 "Use surface layers (1=yes, 0=no, -1=let the algorithm decide). By default (-1), aquamaps will try to find the species' depth range in its internal database. If the minimum depth is equals or less than 200m, then aquamaps will use sea surface layers for temperature and salinity. Otherwise it will use bottom layers. This parameter can be used to force aquamaps to use surface or bottom layers.", // Description 00070 1, // Not zero if the parameter has lower limit 00071 -1, // Parameter's lower limit 00072 1, // Not zero if the parameter has upper limit 00073 1, // Parameter's upper limit 00074 "-1" // Parameter's typical (default) value 00075 }, 00076 { 00077 PARAM_USE_DEPTH_RANGE, // Id 00078 "Use depth range", // Name 00079 Integer, // Type 00080 "Use depth range when calculating probabilities", // Overview 00081 "Use depth range provided by experts (if available) when calculating probabilities", // Description 00082 1, // Not zero if the parameter has lower limit 00083 0, // Parameter's lower limit 00084 1, // Not zero if the parameter has upper limit 00085 1, // Parameter's upper limit 00086 "1" // Parameter's typical (default) value 00087 }, 00088 { 00089 PARAM_USE_ICE_CONCENTRATION, // Id 00090 "Use ice concentration", // Name 00091 Integer, // Type 00092 "Use ice concentration envelope when calculating probabilities", // Overview 00093 "Use ice concentration when calculating probabilities", // Description 00094 1, // Not zero if the parameter has lower limit 00095 0, // Parameter's lower limit 00096 1, // Not zero if the parameter has upper limit 00097 1, // Parameter's upper limit 00098 "1" // Parameter's typical (default) value 00099 }, 00100 { 00101 PARAM_USE_DISTANCE_TO_LAND, // Id 00102 "Use distance to land", // Name 00103 Integer, // Type 00104 "Use distance to land envelope when calculating probabilities", // Overview 00105 "Use distance to land envelope when calculating probabilities", // Description 00106 1, // Not zero if the parameter has lower limit 00107 0, // Parameter's lower limit 00108 1, // Not zero if the parameter has upper limit 00109 1, // Parameter's upper limit 00110 "0" // Parameter's typical (default) value 00111 }, 00112 { 00113 PARAM_USE_PRIMARY_PRODUCTION, // Id 00114 "Use primary production", // Name 00115 Integer, // Type 00116 "Use primary production envelope when calculating probabilities", // Overview 00117 "Use primary production envelope when calculating probabilities", // Description 00118 1, // Not zero if the parameter has lower limit 00119 0, // Parameter's lower limit 00120 1, // Not zero if the parameter has upper limit 00121 1, // Parameter's upper limit 00122 "1" // Parameter's typical (default) value 00123 }, 00124 { 00125 PARAM_USE_SALINITY, // Id 00126 "Use salinity", // Name 00127 Integer, // Type 00128 "Use salinity envelope when calculating probabilities", // Overview 00129 "Use salinity envelope when calculating probabilities", // Description 00130 1, // Not zero if the parameter has lower limit 00131 0, // Parameter's lower limit 00132 1, // Not zero if the parameter has upper limit 00133 1, // Parameter's upper limit 00134 "1" // Parameter's typical (default) value 00135 }, 00136 { 00137 PARAM_USE_TEMPERATURE, // Id 00138 "Use temperature", // Name 00139 Integer, // Type 00140 "Use temperature envelope when calculating probabilities", // Overview 00141 "Use temperature envelope when calculating probabilities", // Description 00142 1, // Not zero if the parameter has lower limit 00143 0, // Parameter's lower limit 00144 1, // Not zero if the parameter has upper limit 00145 1, // Parameter's upper limit 00146 "1" // Parameter's typical (default) value 00147 }, 00148 }; 00149 00150 00151 /************************************/ 00152 /*** Algorithm's general metadata ***/ 00153 00154 static AlgMetadata metadata = { 00155 00156 "AQUAMAPS", // Id. 00157 "AquaMaps (beta version)", // Name. 00158 "0.3", // Version. 00159 00160 // Overview 00161 "Environmental envelope modelling algorithm for marine organisms. \ 00162 Needs 9 predefined layers (see detailed description) to calculate the \ 00163 preferred and accepted envelopes, and also needs a scientific name (genus \ 00164 plus species) labeling each occurrence group." , 00165 00166 // Description. 00167 "AquaMaps is a niche modelling tool developed as part of the Incofish \ 00168 project and particularly tailored towards modelling marine organisms distribution.\n \ 00169 The basic idea follows the environmental envelope type modelling approach, \ 00170 where each variable has an associated preferred range and a broader accepted \ 00171 range. Within the preferred range the probabilty of presence is 1, between the \ 00172 preferred range and the acceptable range the probability varies from 1 to 0 \ 00173 (linear decay), and outside the accepted range the probability is 0. The \ 00174 overall probability is calculated by multiplying all individual probabilities.\n\ 00175 This algorithm differs from other traditional ones since it requires a specific \ 00176 set of layers to work, which should also be in this order: maximum depth in meters, \ 00177 minimum depth in meters, mean annual sea ice concentration, mean annual distance to \ 00178 land in kilometers, mean annual primary production (chlorophyll A, measured in \ 00179 mgC per square meter per day), mean annual bottom salinity in psu, mean annual \ 00180 surface salinity in psu, mean annual bottom temperature in celsius, mean annual \ 00181 surface temperature in celsius. These layers can be downloaded from: \n\ 00182 http://openmodeller.cria.org.br/download/marine2.zip \n\ 00183 Preferred ranges are usually calculated based on percentiles 10th and 90th. \ 00184 They are further adjusted using interquartile values but also ensuring a minimum \ 00185 envelope size based on pre-defined values. Under certain circumstances, AquaMaps can \ 00186 make use of expert information to define the envelopes. This can happen for mammals \ 00187 (any envelope) or for fishes (depth range envelope). Expert information comes from \ 00188 FishBase and is stored in a local SQLite database (aquamaps.db) accessed by the \ 00189 algorithm. To find information in the database, all occurrences provided to \ 00190 openModeller must be labeled with the scientific name (only genus and species). \ 00191 Matches are exact, using case-sensitive equals operation. In this version, the \ 00192 internal database contains information for more then 7000 marine species.", 00193 00194 "Kaschner, K., J. S. Ready, E. Agbayani, J. Rius, K. Kesner-Reyes, P. D. Eastwood, A. B. South, S. O. Kullander, T. Rees, C. H. Close, R. Watson, D. Pauly, and R. Froese", // Authors 00195 00196 // Bibliography. 00197 "Kaschner, K., J. S. Ready, E. Agbayani, J. Rius, K. Kesner-Reyes, P. D. Eastwood, A. B. South, S. O. Kullander, T. Rees, C. H. Close, R. Watson, D. Pauly, and R. Froese. 2007 AquaMaps: Predicted range maps for aquatic species. World wide web electronic publication, www.aquamaps.org, Version 12/2007.", 00198 00199 "Renato De Giovanni", // Code author. 00200 "renato [at] cria dot org dot br", // Code author's contact. 00201 00202 0, // Does not accept categorical data. 00203 0, // Does not need (pseudo)absence points. 00204 00205 NUM_PARAM, // Algorithm's parameters. 00206 parameters 00207 }; 00208 00209 00210 00211 /****************************************************************/ 00212 /****************** Algorithm's factory function ****************/ 00213 00214 OM_ALG_DLL_EXPORT 00215 AlgorithmImpl * 00216 algorithmFactory() 00217 { 00218 return new AquaMaps(); 00219 } 00220 00221 OM_ALG_DLL_EXPORT 00222 AlgMetadata const * 00223 algorithmMetadata() 00224 { 00225 return &metadata; 00226 } 00227 00228 00229 /****************************************************************/ 00230 /**************************** AquaMaps **************************/ 00231 00232 /*******************/ 00233 /*** constructor ***/ 00234 00235 AquaMaps::AquaMaps() : 00236 AlgorithmImpl( &metadata ), 00237 _use_layer( NULL ), 00238 _lower_limit( 7, SURFACE_LOWER_LIMIT ), 00239 _upper_limit( 7, SURFACE_UPPER_LIMIT ), 00240 _inner_size( 7, INNER_SIZE ), 00241 _minimum(), 00242 _maximum(), 00243 _pref_minimum(), 00244 _pref_maximum(), 00245 _pelagic(-1), 00246 _use_surface_layers(-1), 00247 _has_expert_range( NULL ), 00248 _progress( 0.0 ) 00249 { 00250 } 00251 00252 /******************/ 00253 /*** destructor ***/ 00254 00255 AquaMaps::~AquaMaps() 00256 { 00257 if ( _use_layer ) { 00258 00259 delete[] _use_layer; 00260 } 00261 00262 if ( _has_expert_range ) { 00263 00264 delete[] _has_expert_range; 00265 } 00266 } 00267 00268 /******************/ 00269 /*** initialize ***/ 00270 int 00271 AquaMaps::initialize() 00272 { 00273 _has_expert_range = new bool[7]; 00274 00275 // Initialize array 00276 for ( int i = 0; i < 7; ++i ) { 00277 00278 _has_expert_range[i] = false; 00279 } 00280 00281 // Parameter UseSurfaceLayers 00282 if ( ! getParameter( PARAM_USE_SURFACE_LAYERS, &_use_surface_layers ) ) { 00283 00284 _use_surface_layers = -1; 00285 } 00286 00287 _use_layer = new int[7]; 00288 00289 // Parameter UseDepthRange 00290 if ( ! _getAndCheckParameter( PARAM_USE_DEPTH_RANGE, &_use_layer[MAXDEPTH] ) ) { 00291 00292 return 0; 00293 } 00294 else { 00295 00296 _use_layer[MINDEPTH] = _use_layer[MAXDEPTH]; // just repeat the value 00297 } 00298 00299 // Parameter UseIceConcentration 00300 if ( ! _getAndCheckParameter( PARAM_USE_ICE_CONCENTRATION, &_use_layer[ICE_CONCENTRATION] ) ) { 00301 00302 return 0; 00303 } 00304 00305 // Parameter UseDistanceToLand 00306 if ( ! _getAndCheckParameter( PARAM_USE_DISTANCE_TO_LAND, &_use_layer[DISTANCE_TO_LAND] ) ) { 00307 00308 return 0; 00309 } 00310 00311 // Parameter UsePrimaryProduction 00312 if ( ! _getAndCheckParameter( PARAM_USE_PRIMARY_PRODUCTION, &_use_layer[PRIMARY_PRODUCTION] ) ) { 00313 00314 return 0; 00315 } 00316 00317 // Parameter UseSalinity 00318 if ( ! _getAndCheckParameter( PARAM_USE_SALINITY, &_use_layer[SALINITY] ) ) { 00319 00320 return 0; 00321 } 00322 00323 // Parameter UseTemperature 00324 if ( ! _getAndCheckParameter( PARAM_USE_TEMPERATURE, &_use_layer[TEMPERATURE] ) ) { 00325 00326 return 0; 00327 } 00328 00329 // Number of independent variables. 00330 int dim = _samp->numIndependent(); 00331 Log::instance()->info( "Reading %d-dimensional occurrence points.\n", dim ); 00332 00333 // Check the number of layers. 00334 if ( dim != 9 ) { 00335 00336 Log::instance()->error( "AquaMaps needs precisely 9 layers to work, and they should be in this order: maximum depth in meters, minimum depth in meters, mean annual sea ice concentration, mean annual distance to land in Kilometers, mean annual primary production (chlorophyll A), mean annual sea bottom salinity in psu,, mean annual sea surface salinity in psu, mean annual sea bottom temperature in Celsius, mean annual sea surface temperature in Celsius.\n" ); 00337 return 0; 00338 } 00339 00340 // Remove duplicates accross geography 00341 _samp->spatiallyUnique(); 00342 00343 // Check the number of sampled points. 00344 int npnt = _samp->numPresence(); 00345 00346 if ( npnt < 5 ) { 00347 00348 Log::instance()->error( "AquaMaps needs at least 5 points inside the mask!\n" ); 00349 return 0; 00350 } 00351 00352 return 1; 00353 } 00354 00355 00356 /***************/ 00357 /*** iterate ***/ 00358 int 00359 AquaMaps::iterate() 00360 { 00361 Log::instance()->info( "Using %d points to find AquaMaps envelopes.\n", _samp->numPresence() ); 00362 00363 _calculateEnvelopes( _samp->getPresences() ); 00364 00365 return 1; 00366 } 00367 00368 00369 /***********************/ 00370 /*** get Convergence ***/ 00371 int 00372 AquaMaps::getConvergence( Scalar * const val ) const 00373 { 00374 *val = 1.0; 00375 return 1; 00376 } 00377 00378 00379 /********************/ 00380 /*** get Progress ***/ 00381 float 00382 AquaMaps::getProgress() const 00383 { 00384 return _progress; 00385 } 00386 00387 00388 /************/ 00389 /*** done ***/ 00390 int 00391 AquaMaps::done() const 00392 { 00393 // This is not an iterative algorithm. 00394 return ( _progress > 0.99 ) ? 1 : 0; 00395 } 00396 00397 00398 /*******************************/ 00399 /*** get and check parameter ***/ 00400 int 00401 AquaMaps::_getAndCheckParameter( std::string const &name, int * value ) 00402 { 00403 // Parameter UseSalinity 00404 if ( ! getParameter( name, value ) ) { 00405 00406 Log::instance()->error( "Parameter '%s' not passed.\n", name.c_str() ); 00407 return 0; 00408 } 00409 // Check if value is 0 or 1 00410 if ( *value != 0 && *value != 1 ) { 00411 00412 Log::instance()->error( "Parameter '%s' not set properly. It must be 0 or 1.\n", name.c_str() ); 00413 return 0; 00414 } 00415 00416 return 1; 00417 } 00418 00419 00420 /***************************/ 00421 /*** calculate envelopes ***/ 00422 void 00423 AquaMaps::_calculateEnvelopes( const OccurrencesPtr& occs ) 00424 { 00425 Log::instance()->debug("Species is: %s\n", occs->label()); 00426 Log::instance()->debug("Layers are:\n"); 00427 Log::instance()->debug("0 = Maximum depth\n"); 00428 Log::instance()->debug("1 = Minimum depth\n"); 00429 Log::instance()->debug("2 = Ice concentration\n"); 00430 Log::instance()->debug("3 = Distance to land\n"); 00431 Log::instance()->debug("4 = Primary production (chlorophyll A)\n"); 00432 Log::instance()->debug("5 = Salinity bottom\n"); 00433 Log::instance()->debug("6 = Salinity surface\n"); 00434 Log::instance()->debug("7 = Temperature bottom\n"); 00435 Log::instance()->debug("8 = Temperature surface\n"); 00436 00437 // Compute min, pref_min, pref_max and max 00438 OccurrencesImpl::const_iterator oc = occs->begin(); 00439 OccurrencesImpl::const_iterator ocEnd = occs->end(); 00440 00441 // Intialize _minimum, _maximum, _pref_minimum, and _pref_maximum 00442 // to the values of the first point, and increment 00443 // to get it out of the next loop. 00444 Sample const & sample = (*oc)->environment(); 00445 _minimum = sample; 00446 _pref_minimum = sample; 00447 _pref_maximum = sample; 00448 _maximum = sample; 00449 Sample mean = sample; 00450 ++oc; 00451 00452 // Load default min/max 00453 while ( oc != ocEnd ) { 00454 00455 Sample const& sample = (*oc)->environment(); 00456 00457 // Default min/max will be the actual min/max values 00458 _minimum &= sample; 00459 _maximum |= sample; 00460 00461 mean += sample; 00462 00463 ++oc; 00464 } 00465 00466 mean /= occs->numOccurrences(); 00467 00468 int num_layers = _minimum.size(); 00469 00470 _progress = 1.0f/num_layers; // this is arbitrary 00471 00472 // Default values for depth ranges 00473 _minimum[MAXDEPTH] = _minimum[MINDEPTH] = _pref_minimum[MAXDEPTH] = _pref_minimum[MINDEPTH] = 0.0; 00474 _maximum[MAXDEPTH] = _maximum[MINDEPTH] = _pref_maximum[MAXDEPTH] = _pref_maximum[MINDEPTH] = 9999.0; 00475 00476 // Try to get expert information about depth range from database 00477 _readSpeciesData( occs->label() ); 00478 00479 if ( _use_surface_layers == -1 ) { 00480 00481 _use_surface_layers = ( _minimum[MINDEPTH] <= DEPTH_LIMIT ) ? 0 : 1; 00482 } 00483 00484 if ( ! _use_surface_layers ) { 00485 00486 Log::instance()->info("Using bottom layers.\n"); 00487 00488 for ( int i = 0; i < 7; ++i ) { 00489 00490 _lower_limit[i] = BOTTOM_LOWER_LIMIT[i]; 00491 _upper_limit[i] = BOTTOM_UPPER_LIMIT[i]; 00492 } 00493 } 00494 else { 00495 00496 Log::instance()->info("Using surface layers.\n"); 00497 } 00498 00499 _progress = 2.0f/num_layers; // this is arbitrary 00500 00501 // Get matrix data structure so that we can sort values for each layer 00502 std::vector<ScalarVector> matrix = occs->getEnvironmentMatrix(); 00503 00504 Scalar nodata = -1.0; 00505 00506 // For each layer (except min/max depth), calculate the percentiles and set 00507 // default values for the prefered min/max 00508 ScalarVector::iterator lStart, lEnd; 00509 00510 for ( int j = 2; j < num_layers; j++ ) { 00511 00512 if ( _use_surface_layers ) { 00513 00514 // Ignore bottom layers 00515 if ( j == 5 || j == 7 ) { 00516 00517 _progress = (float)(j+1)/num_layers; 00518 00519 _minimum[j] = _minimum[j] = _pref_minimum[j] = _pref_minimum[j] = nodata; 00520 00521 continue; 00522 } 00523 } 00524 else { 00525 00526 // Ignore surface layers 00527 if ( j == 6 || j == 8 ) { 00528 00529 _progress = (float)(j+1)/num_layers; 00530 00531 _minimum[j] = _minimum[j] = _pref_minimum[j] = _pref_minimum[j] = nodata; 00532 00533 continue; 00534 } 00535 } 00536 00537 if ( ( j == 5 || j == 6 ) && _has_expert_range[SALINITY] ) { 00538 00539 _progress = (float)(j+1)/num_layers; 00540 00541 continue; 00542 } 00543 00544 if ( ( j == 7 || j == 8 ) && _has_expert_range[TEMPERATURE] ) { 00545 00546 _progress = (float)(j+1)/num_layers; 00547 00548 continue; 00549 } 00550 00551 Log::instance()->debug("--------------------------------\n", j); 00552 Log::instance()->debug("Calculating envelope for layer %d\n", j); 00553 Log::instance()->debug("--------------------------------\n", j); 00554 00555 // 0 = Maximum depth 00556 // 1 = Minimum depth 00557 // 2 = Ice concentration 00558 // 3 = Distance to land 00559 // 4 = Primary production (chlorophyll A) 00560 // 5 = Salinity bottom 00561 // 6 = Salinity surface 00562 // 7 = Temperature bottom 00563 // 8 = Temperature surface 00564 00565 lStart = matrix[j].begin(); 00566 lEnd = matrix[j].end(); 00567 sort( lStart, lEnd ); 00568 00569 int numOccurrences = occs->numOccurrences(); 00570 00571 // Just for debugging 00572 //ostringstream debug; 00573 //debug << "row -->"; 00574 //for ( int w = 0; w < numOccurrences; w++ ) { 00575 00576 //debug << " [" << j << "," << w << "]=" << matrix[j][w]; 00577 //} 00578 //debug << "\n"; 00579 //Log::instance()->debug( debug.str().c_str() ); 00580 00581 // Calculate percentiles 00582 Scalar v10, v25, v75, v90; 00583 00584 _percentile( &v10, numOccurrences, 0.10, &matrix, j ); 00585 _percentile( &v25, numOccurrences, 0.25, &matrix, j ); 00586 _percentile( &v75, numOccurrences, 0.75, &matrix, j ); 00587 _percentile( &v90, numOccurrences, 0.90, &matrix, j ); 00588 00589 Log::instance()->debug("10th percentile: %f\n", v10); 00590 Log::instance()->debug("25th percentile: %f\n", v25); 00591 Log::instance()->debug("75th percentile: %f\n", v75); 00592 Log::instance()->debug("90th percentile: %f\n", v90); 00593 00594 Scalar interquartile = fabs( v25 - v75 ); 00595 00596 Log::instance()->debug("Interquartile: %f\n", interquartile); 00597 00598 Scalar adjmin = v25 - ( 1.5 * interquartile ); 00599 Scalar adjmax = v75 + ( 1.5 * interquartile ); 00600 00601 Log::instance()->debug("Adjmin: %f\n", adjmin); 00602 Log::instance()->debug("Adjmax: %f\n", adjmax); 00603 00604 _pref_minimum[j] = v10; 00605 _pref_maximum[j] = v90; 00606 00607 Log::instance()->debug("_Before adjustments_\n"); 00608 Log::instance()->debug("min: %f\n", _minimum[j]); 00609 Log::instance()->debug("prefmin: %f\n", _pref_minimum[j]); 00610 Log::instance()->debug("prefmax: %f\n", _pref_maximum[j]); 00611 Log::instance()->debug("max: %f\n", _maximum[j]); 00612 00613 // Make the interquartile adjusting and ensure the envelope sizes 00614 // for all variables except ice concentration 00615 if ( j > 2 ) { 00616 00617 _adjustInterquartile( j, adjmin, adjmax ); 00618 _ensureEnvelopeSize( j ); 00619 } 00620 else { 00621 00622 if ( _minimum[j] == 0.0) { 00623 00624 _minimum[j] = mean[j] - 1.0; // black magic copied from original vb code 00625 } 00626 } 00627 00628 Log::instance()->debug("_After adjustments_\n"); 00629 Log::instance()->debug("min: %f\n", _minimum[j]); 00630 Log::instance()->debug("prefmin: %f\n", _pref_minimum[j]); 00631 Log::instance()->debug("prefmax: %f\n", _pref_maximum[j]); 00632 Log::instance()->debug("max: %f\n", _maximum[j]); 00633 00634 _progress = (float)(j+1)/num_layers; 00635 } 00636 } 00637 00638 00639 /******************/ 00640 /*** percentile ***/ 00641 void 00642 AquaMaps::_percentile( Scalar *result, int n, double percent, std::vector<ScalarVector> *matrix, int layerIndex ) 00643 { 00644 int iPrev, iNext; 00645 00646 Scalar k = n * percent; 00647 00648 iPrev = (int) floor( k ); 00649 iNext = (int) ceil( k ); 00650 00651 // Subtract 1 from the result since arrays begin with position 0 00652 --iPrev; 00653 --iNext; 00654 00655 if ( iPrev == iNext ) { 00656 00657 // If k is an integer, use the mean between the two positions 00658 *result = ( (*matrix)[layerIndex][iPrev] + (*matrix)[layerIndex][iNext] ) / 2; 00659 } 00660 else { 00661 00662 // If k is not an integer, use the next as a reference 00663 *result = (*matrix)[layerIndex][iNext]; 00664 } 00665 } 00666 00667 00668 /***********************/ 00669 /*** read depth data ***/ 00670 void 00671 AquaMaps::_readSpeciesData( const char *species ) 00672 { 00673 string dbfullname = omDataPath(); 00674 00675 #ifndef WIN32 00676 dbfullname.append( "/" ); 00677 #else 00678 dbfullname.append( "\\" ); 00679 #endif 00680 00681 dbfullname.append( "aquamaps.db" ); 00682 00683 Log::instance()->debug( "Trying database: %s\n", dbfullname.c_str() ); 00684 00685 ifstream dbfile( dbfullname.c_str(), ios::in ); 00686 00687 if ( ! dbfile ) { 00688 00689 Log::instance()->warn( "Could not find internal database with species data.\n" ); 00690 return; 00691 } 00692 00693 Log::instance()->info( "Using internal database with species data.\n" ); 00694 00695 sqlite3 *db; 00696 00697 int rc = sqlite3_open( dbfullname.c_str(), &db); 00698 00699 if ( rc ) { 00700 00701 // This will likely never happen since on open, sqlite creates the 00702 // database if it does not exist. 00703 Log::instance()->warn( "Could not open internal database: %s\n", sqlite3_errmsg( db ) ); 00704 sqlite3_close( db ); 00705 return; 00706 } 00707 00708 // Prepare the sql statement 00709 const char *pzTail; 00710 00711 sqlite3_stmt *ppStmt; 00712 00713 char* sql = sqlite3_mprintf( "select pelagic, provider, depthmin, depthmax, depthprefmin, depthprefmax, iceconmin, iceconmax, iceconprefmin, iceconprefmax, landdistmin, landdistmax, landdistprefmin, landdistprefmax, primprodmin, primprodmax, primprodprefmin, primprodprefmax, salinitymin, salinitymax, salinityprefmin, salinityprefmax, tempmin, tempmax, tempprefmin, tempprefmax from spinfo where species = '%q'", species ); 00714 00715 rc = sqlite3_prepare( db, sql, strlen( sql ), &ppStmt, &pzTail ); 00716 00717 if ( rc == SQLITE_OK ) { 00718 00719 // Fecth data 00720 int result_code = sqlite3_step( ppStmt ); 00721 00722 if ( result_code == SQLITE_ROW ) { 00723 00724 int pelagic = sqlite3_column_int( ppStmt, 0 ); 00725 const char * provider = (char *)sqlite3_column_text( ppStmt, 1 ); 00726 double depthmin = sqlite3_column_double( ppStmt, 2 ); 00727 double depthmax = sqlite3_column_double( ppStmt, 3 ); 00728 double depthprefmin = sqlite3_column_double( ppStmt, 4 ); 00729 double depthprefmax = sqlite3_column_double( ppStmt, 5 ); 00730 00731 _pelagic = pelagic; 00732 00733 // Repeat the information for both min and max depth for symmetry 00734 _minimum[MAXDEPTH] = depthmin; 00735 _minimum[MINDEPTH] = depthmin; 00736 _maximum[MAXDEPTH] = depthmax; 00737 _maximum[MINDEPTH] = depthmax; 00738 _pref_minimum[MAXDEPTH] = depthprefmin; 00739 _pref_minimum[MINDEPTH] = depthprefmin; 00740 _pref_maximum[MAXDEPTH] = depthprefmax; 00741 _pref_maximum[MINDEPTH] = depthprefmax; 00742 00743 Log::instance()->debug("Values from expert database:\n"); 00744 Log::instance()->debug("provider: %s\n", provider); 00745 Log::instance()->debug("pelagic: %i\n", pelagic); 00746 Log::instance()->debug("depthmin: %f\n", depthmin); 00747 Log::instance()->debug("depthmax: %f\n", depthmax); 00748 Log::instance()->debug("depthprefmin: %f\n", depthprefmin); 00749 Log::instance()->debug("depthprefmax: %f\n", depthprefmax); 00750 00751 if ( sqlite3_column_bytes( ppStmt, 1 ) ) { 00752 00753 string prov_code( provider ); 00754 00755 if ( prov_code == "MM" ) { 00756 00757 Log::instance()->info( "Looks like a mammal\n" ); 00758 00759 for ( int i = ICE_CONCENTRATION; i < TEMPERATURE + 1; ++i ) { 00760 00761 _has_expert_range[i] = _hasExpertRange( ppStmt, i ); 00762 00763 if ( _has_expert_range[i] ) { 00764 00765 Log::instance()->info( "Found expert range for variable %s\n", NAME[i].c_str() ); 00766 00767 Log::instance()->debug("min: %f\n", _minimum[i]); 00768 Log::instance()->debug("max: %f\n", _maximum[i]); 00769 Log::instance()->debug("prefmin: %f\n", _pref_minimum[i]); 00770 Log::instance()->debug("prefmax: %f\n", _pref_maximum[i]); 00771 } 00772 } 00773 } 00774 } 00775 } 00776 else if ( result_code == SQLITE_DONE ) { 00777 00778 Log::instance()->warn( "'%s' not found in species database\n", species ); 00779 } 00780 else { 00781 00782 Log::instance()->warn( "Could not fetch data from species database: %s\n", sqlite3_errmsg( db ) ); 00783 } 00784 } 00785 else { 00786 00787 Log::instance()->warn( "Could not prepare SQL statement to query species database: %s\n", sqlite3_errmsg( db ) ); 00788 } 00789 00790 // free memory 00791 sqlite3_free( sql ); 00792 00793 // close the statement 00794 sqlite3_finalize( ppStmt ); 00795 00796 // close the database 00797 sqlite3_close( db ); 00798 } 00799 00800 00801 /****************************/ 00802 /*** ensure envelope size ***/ 00803 bool 00804 AquaMaps::_hasExpertRange( sqlite3_stmt * stmt, int varIndex ) 00805 { 00806 int index; 00807 00808 if ( varIndex == ICE_CONCENTRATION ) { 00809 00810 index = 6; 00811 } 00812 else if ( varIndex == DISTANCE_TO_LAND ) { 00813 00814 index = 10; 00815 } 00816 else if ( varIndex == PRIMARY_PRODUCTION ) { 00817 00818 index = 14; 00819 } 00820 else if ( varIndex == SALINITY ) { 00821 00822 index = 18; 00823 } 00824 else if ( varIndex == TEMPERATURE ) { 00825 00826 index = 22; 00827 } 00828 else { 00829 00830 return false; 00831 } 00832 00833 if ( ! sqlite3_column_bytes( stmt, index ) ) { 00834 00835 return false; 00836 } 00837 00838 double value = sqlite3_column_double( stmt, index ); 00839 00840 _minimum[varIndex] = value; 00841 00842 ++index; 00843 00844 if ( ! sqlite3_column_bytes( stmt, index ) ) { 00845 00846 return false; 00847 } 00848 00849 value = sqlite3_column_double( stmt, index ); 00850 00851 _maximum[varIndex] = value; 00852 00853 ++index; 00854 00855 if ( ! sqlite3_column_bytes( stmt, index ) ) { 00856 00857 return false; 00858 } 00859 00860 value = sqlite3_column_double( stmt, index ); 00861 00862 _pref_minimum[varIndex] = value; 00863 00864 ++index; 00865 00866 if ( ! sqlite3_column_bytes( stmt, index ) ) { 00867 00868 return false; 00869 } 00870 00871 value = sqlite3_column_double( stmt, index ); 00872 00873 _pref_maximum[varIndex] = value; 00874 00875 return true; 00876 } 00877 00878 00879 /****************************/ 00880 /*** ensure envelope size ***/ 00881 void 00882 AquaMaps::_adjustInterquartile( int layerIndex, Scalar adjmin, Scalar adjmax ) 00883 { 00884 if ( adjmin > _lower_limit[_getRelatedIndex(layerIndex)] && adjmin < _minimum[layerIndex] ) { 00885 00886 _minimum[layerIndex] = adjmin; 00887 } 00888 00889 if ( adjmax < _upper_limit[_getRelatedIndex(layerIndex)] && adjmax > _maximum[layerIndex] ) { 00890 00891 _maximum[layerIndex] = adjmax; 00892 } 00893 } 00894 00895 00896 /****************************/ 00897 /*** ensure envelope size ***/ 00898 void 00899 AquaMaps::_ensureEnvelopeSize( int layerIndex ) 00900 { 00901 Scalar halfSize = _inner_size[_getRelatedIndex(layerIndex)]/2; 00902 00903 // Try to ensure that the prefered envelope has the specified size 00904 // (but only change if the new values are still within the accepted envelope) 00905 if ( _pref_maximum[layerIndex] - _pref_minimum[layerIndex] < _inner_size[_getRelatedIndex(layerIndex)] ) { 00906 00907 Scalar center = (_pref_minimum[layerIndex] + _pref_maximum[layerIndex]) / 2; 00908 00909 Scalar new_pref_min = center - halfSize; 00910 00911 if ( new_pref_min >= _minimum[layerIndex] ) { 00912 00913 _pref_minimum[layerIndex] = new_pref_min; 00914 } 00915 00916 Scalar new_pref_max = center + halfSize; 00917 00918 if ( new_pref_max <= _maximum[layerIndex] ) { 00919 00920 _pref_maximum[layerIndex] = new_pref_max; 00921 } 00922 } 00923 00924 // Ensure that the distance between the envelopes has half the specified size 00925 // (but changes should never exceed the allowed limits!) 00926 if ( _pref_minimum[layerIndex] - _minimum[layerIndex] < halfSize ) { 00927 00928 Scalar new_min = _pref_minimum[layerIndex] - halfSize; 00929 00930 if ( new_min > _lower_limit[_getRelatedIndex(layerIndex)] ) { 00931 00932 _minimum[layerIndex] = new_min; 00933 } 00934 else { 00935 00936 _minimum[layerIndex] = _lower_limit[_getRelatedIndex(layerIndex)]; 00937 } 00938 } 00939 00940 if ( _maximum[layerIndex] - _pref_maximum[layerIndex] < halfSize ) { 00941 00942 Scalar new_max = _pref_maximum[layerIndex] + halfSize; 00943 00944 if ( new_max < _upper_limit[_getRelatedIndex(layerIndex)] ) { 00945 00946 _maximum[layerIndex] = new_max; 00947 } 00948 else { 00949 00950 _maximum[layerIndex] = _upper_limit[_getRelatedIndex(layerIndex)]; 00951 } 00952 } 00953 } 00954 00955 00956 /*****************/ 00957 /*** get Value ***/ 00958 Scalar 00959 AquaMaps::getValue( const Sample& x ) const 00960 { 00961 // Final probability is the multiplication of all 00962 // individual probabilities 00963 Scalar prob = 1.0; 00964 00965 int num_variables_used = 0; 00966 00967 // Depth probability 00968 00969 if ( _use_layer[MAXDEPTH] ) { // If user wants to use depth range 00970 00971 ++num_variables_used; 00972 00973 if ( _maximum[MAXDEPTH] != 9999 ) { // If there is a depth range in the internal database 00974 00975 // Probability of occurrence is zero if depth at this point is less 00976 // than the minimum depth for the species. 00977 // note: using maximum depth layer here 00978 if ( x[MAXDEPTH] < _minimum[MAXDEPTH] ) { 00979 00980 return 0.0; 00981 } 00982 00983 // Point is sufficiently deep, but still not enough for the prefered range 00984 // note: using maximum depth layer here 00985 if ( x[MAXDEPTH] >= _minimum[MAXDEPTH] && x[MAXDEPTH] < _pref_minimum[MAXDEPTH] ) { 00986 00987 // Linear decay 00988 prob *= (x[MAXDEPTH] - _minimum[MAXDEPTH]) / (_pref_minimum[MAXDEPTH] - _minimum[MAXDEPTH]); 00989 } 00990 // Point is sufficiently deep for the prefered range 00991 // Note: pelagic means "belonging to the upper layers of the open sea". Or in other 00992 // words, no need to feed at the bottom of the sea. 00993 else if ( _pelagic == 1 ) { 00994 00995 // Just keep prob as 1 00996 } 00997 // Species needs to feed at the bottom of the sea (or there's no data about this) 00998 // and point is within the prefered range 00999 // note: the inner envelope logic makes use of both minimum and maximum depth layers. 01000 else if ( x[MAXDEPTH] >= _pref_minimum[MAXDEPTH] && x[MINDEPTH] <= _pref_maximum[MINDEPTH] ) { 01001 01002 // Just keep prob as 1 01003 } 01004 // Species needs to feed at the bottom of the sea (or there's no data about this) and 01005 // point is deeper than the prefered maximum depth but not so deep as the maximum 01006 // note: using minimum depth layer here 01007 else if ( x[MINDEPTH] >= _pref_maximum[MINDEPTH] && x[MINDEPTH] < _maximum[MINDEPTH] ) { 01008 01009 // Linear decay 01010 prob *= (_maximum[MINDEPTH] - x[MINDEPTH]) / (_maximum[MINDEPTH] - _pref_maximum[MINDEPTH]); 01011 } 01012 // Species needs to feed at the bottom of the sea (or there's no data about this) 01013 // but depth at the point is greater then the maximum accepted depth 01014 else { 01015 01016 return 0.0; 01017 } 01018 } 01019 else { // If there is no depth range 01020 01021 // Just keep prob as 1 01022 } 01023 } 01024 else { // If user doesn't want to use depth ranges 01025 01026 // Just keep prob as 1 01027 } 01028 01029 // For all layers except min/max depths 01030 01031 int numLayers = x.size(); 01032 01033 for ( int i = 2; i < numLayers; i++ ) { 01034 01035 // Ignore layer if user doesn't want to use it 01036 if ( i < 5 ) { 01037 01038 if ( ! _use_layer[i] ) { 01039 01040 continue; 01041 } 01042 } 01043 else if ( i < 7 ) { 01044 01045 if ( ! _use_layer[SALINITY] ) { 01046 01047 continue; 01048 } 01049 } 01050 else { 01051 01052 if ( ! _use_layer[TEMPERATURE] ) { 01053 01054 continue; 01055 } 01056 } 01057 01058 if ( _use_surface_layers ) { 01059 01060 // Ignore bottom layers 01061 if ( i == 5 || i == 7 ) { 01062 01063 continue; 01064 } 01065 } 01066 else 01067 { 01068 // Ignore surface layers 01069 if ( i == 6 || i == 8 ) { 01070 01071 continue; 01072 } 01073 } 01074 01075 ++num_variables_used; 01076 01077 // Probability zero for points outside the envelope 01078 if ( x[i] < _minimum[i] || x[i] > _maximum[i] ) { 01079 01080 return 0.0; 01081 } 01082 01083 // Linear decay for points outside the prefered envelope 01084 if ( x[i] < _pref_minimum[i] ) { 01085 01086 prob *= (x[i] - _minimum[i]) / (_pref_minimum[i] - _minimum[i]); 01087 } 01088 // Linear decay for points outside the prefered envelope 01089 else if ( x[i] > _pref_maximum[i] ) { 01090 01091 prob *= (_maximum[i] - x[i]) / (_maximum[i] - _pref_maximum[i]); 01092 } 01093 else { 01094 01095 // No need to check inside preferred envelope - probability is 1 there, so 01096 // it would not change the final value 01097 } 01098 } 01099 01100 // Return probability product 01101 //return prob; 01102 01103 // Return geometric mean 01104 return pow( prob, (Scalar)1/num_variables_used ); 01105 } 01106 01107 01108 /****************************************************************/ 01109 /****************** configuration *******************************/ 01110 void 01111 AquaMaps::_getConfiguration( ConfigurationPtr& config ) const 01112 { 01113 if ( ! done() ) { 01114 01115 return; 01116 } 01117 01118 ConfigurationPtr model_config( new ConfigurationImpl( "AquaMaps" ) ); 01119 config->addSubsection( model_config ); 01120 01121 model_config->addNameValue( "UseLayer", _use_layer, 7 ); 01122 model_config->addNameValue( "UseSurfaceLayers", _use_surface_layers ); 01123 model_config->addNameValue( "Pelagic", _pelagic ); 01124 model_config->addNameValue( "Minimum", _minimum ); 01125 model_config->addNameValue( "PreferedMinimum", _pref_minimum ); 01126 model_config->addNameValue( "PreferedMaximum", _pref_maximum ); 01127 model_config->addNameValue( "Maximum", _maximum ); 01128 } 01129 01130 01131 void 01132 AquaMaps::_setConfiguration( const ConstConfigurationPtr& config ) 01133 { 01134 ConstConfigurationPtr model_config = config->getSubsection( "AquaMaps", false ); 01135 01136 if ( ! model_config ) { 01137 01138 return; 01139 } 01140 01141 try { 01142 01143 int size; 01144 01145 model_config->getAttributeAsIntArray( "UseLayer", &_use_layer, &size ); 01146 01147 _use_surface_layers = model_config->getAttributeAsInt( "UseSurfaceLayers", -1 ); 01148 } 01149 catch ( AttributeNotFound& exception ) { 01150 01151 Log::instance()->error( "Could not find attribute '%s' in serialized model. You may be trying to load a model that was generated with an older version (< 0.2) of this algorithm.\n", exception.getName().c_str() ); 01152 01153 return; 01154 } 01155 01156 _pelagic = model_config->getAttributeAsInt( "Pelagic", -1 ); 01157 _minimum = model_config->getAttributeAsSample( "Minimum" ); 01158 _pref_minimum = model_config->getAttributeAsSample( "PreferedMinimum" ); 01159 _pref_maximum = model_config->getAttributeAsSample( "PreferedMaximum" ); 01160 _maximum = model_config->getAttributeAsSample( "Maximum" ); 01161 01162 _progress = 1.0; 01163 01164 return; 01165 } 01166 01167 01168 /********************/ 01169 /*** log Envelope ***/ 01170 void 01171 AquaMaps::_logEnvelope() 01172 { 01173 Log::instance()->info( "Envelope with %d dimensions (variables).\n\n", _minimum.size() ); 01174 01175 for ( unsigned int i = 0; i < _minimum.size(); i++ ) { 01176 01177 Log::instance()->info( "Variable %02d:", i ); 01178 Log::instance()->info( " Minimum : %f\n", _minimum[i] ); 01179 Log::instance()->info( " Prefered Minimum: %f\n", _pref_minimum[i] ); 01180 Log::instance()->info( " Prefered Maximum: %f\n", _pref_maximum[i] ); 01181 Log::instance()->info( " Maximum : %f\n", _maximum[i] ); 01182 Log::instance()->info( "\n" ); 01183 } 01184 } 01185 01186 int 01187 AquaMaps::_getRelatedIndex( int index ) 01188 { 01189 if ( index == 8 ) { 01190 01191 return index - 2; 01192 } 01193 else if ( index > 5 ) { 01194 01195 return index - 1; 01196 } 01197 01198 return index; 01199 }