Example demonstrating continent boundary detection on spherical meshes. More...

#include "moab/Core.hpp"
#include "moab/ReadUtilIface.hpp"
#include <cmath>
#include "moab/CartVect.hpp"
#include <iostream>
#include <fstream>

Include dependency graph for ContinentsOnGlobe.cpp:

Go to the source code of this file.

Functions
double	getLat (CartVect p)

double	getLon (CartVect p)

bool	interior_point (vector< double > &coords, int &startLoop, int &endLoop, double lat, double lon)

int	main (int argc, char **argv)

Variables
string	bd_name = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" )

string	loops = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" )

string	input_file = string( MESH_DIR ) + string( "/../examples/earth/poly2000.h5m" )

Detailed Description

Example demonstrating continent boundary detection on spherical meshes.

This example shows how to:

Load spherical mesh files and continent boundary data
Convert 3D spherical coordinates to 2D lat-lon coordinates
Perform point-in-polygon tests for continent detection
Create mesh sets for different continents and islands
Tag mesh elements with continent information
Handle complex boundary loops for major landmasses
Write continent-mapped mesh files

This tool is useful for climate and geophysical applications where continent boundaries need to be identified on spherical meshes.

Author: MOAB Development Team

Date: 2024

Definition in file ContinentsOnGlobe.cpp.

Function Documentation

◆ getLat()

double getLat ( CartVect p )

Examples: ContinentsOnGlobe.cpp.

Definition at line 49 of file ContinentsOnGlobe.cpp.

 {
     p.normalize();
     return asin( p[2] );
 }

References moab::CartVect::normalize().

Referenced by main().

◆ getLon()

double getLon ( CartVect p )

Examples: ContinentsOnGlobe.cpp.

Definition at line 54 of file ContinentsOnGlobe.cpp.

 {
     p.normalize();
     double lon;
  
     lon = atan2( p[1], p[0] );
  
     if( lon < -2.95 )
     {  // separate Asia/Europe from Americas
         return 2.0 * M_PI + lon;
     }
     else
     {
         return lon;
     }
 }

References moab::CartVect::normalize().

Referenced by main().

◆ interior_point()

bool interior_point	(	vector< double > &	coords,
		int &	startLoop,
		int &	endLoop,
		double	lat,
		double	lon
	)

Examples: ContinentsOnGlobe.cpp.

Definition at line 74 of file ContinentsOnGlobe.cpp.

 {
     // compute oriented angle between point and edges in the loop
     // if angle ~ +/-2Pi, it is in the interior
     double totalAngle = 0;
     for( int i = startLoop; i <= endLoop; i++ )
     {
         double x1 = coords[2 * i], y1 = coords[2 * i + 1];
         int i2 = i + 1;
         if( i == endLoop ) i2 = startLoop;
         double x2 = coords[2 * i2], y2 = coords[2 * i2 + 1];
         CartVect P1( x1 - lat, y1 - lon, 0. );
         CartVect P2( x2 - lat, y2 - lon, 0. );
         CartVect cross = P1 * P2;
         double angle1  = angle( P1, P2 );
         if( cross[2] < 0 ) angle1 = -angle1;
         totalAngle += angle1;
     }
     if( fabs( totalAngle + 2 * M_PI ) < 0.1 || fabs( totalAngle - 2 * M_PI ) < 0.1 ) return true;
     return false;
 }

References moab::angle(), and moab::cross().

Referenced by main().

◆ main()

int main	(	int	argc,
		char **	argv
	)

get loops for boundaries

Examples: ContinentsOnGlobe.cpp.

Definition at line 96 of file ContinentsOnGlobe.cpp.

 {
     // Get MOAB instance
     Interface* mb = new( std::nothrow ) Core;
     if( NULL == mb ) return 1;
  
     if( argc > 1 )
     {
         // User has input file for mesh file
         input_file = argv[1];
     }
  
     cout << "input file : " << input_file << "\n";
     std::vector< double > coords;
     ifstream bdFile( bd_name.c_str() );
     if( !bdFile )
     {
         cout << endl << "Failed to open file " << bd_name << endl;
         ;
         return 1;
     }
     coords.resize( 0 );
     while( !bdFile.eof() )
     {
         double co;
         bdFile >> co;
         coords.push_back( co );
     }
     cout << " number of boundary points:" << coords.size() / 3 << "\n";
     /// get loops for boundaries
     vector< int > loopsindx;
     ifstream loopFile( loops.c_str() );
     while( !loopFile.eof() )
     {
         int indx;
         loopFile >> indx;
         loopsindx.push_back( indx );
     }
     cout << "number of loops: " << loopsindx.size() / 2 << "\n";
     // convert loops to lat-lon loops
     // it will be easier to determine in-out relations for points in 2d plane
     // careful with the international time zone; ( + - Pi longitude )
     vector< double > coords2d;
     coords2d.resize( coords.size() / 3 * 2 );
     for( int i = 0; i < (int)coords.size() / 3; i++ )
     {
         CartVect p( coords[3 * i], coords[3 * i + 1], coords[3 * i + 2] );
         coords2d[2 * i]     = getLat( p );
         coords2d[2 * i + 1] = getLon( p );
     }
  
     MB_CHK_SET_ERR( mb->load_file( input_file.c_str() ), "Can't load file" );
     // look at the center of element, and see if it is inside the loop
  
     Range cells;
     MB_CHK_SET_ERR( mb->get_entities_by_dimension( 0, 2, cells ), "Can't get cells" );
  
     cout << "number of cells: " << cells.size() << "\n";
  
     // tag for continents
     Tag tag1;
     int defa = -1;
     MB_CHK_SET_ERR( mb->tag_get_handle( "continent", 1, MB_TYPE_INTEGER, tag1, MB_TAG_DENSE | MB_TAG_CREAT, &defa ),
                     "Trouble creating continent tag" );
     EntityHandle islandSets[6];
     for( int loop_index = 0; loop_index < 6; loop_index++ )
     {
         MB_CHK_SET_ERR( mb->create_meshset( MESHSET_SET, islandSets[loop_index] ), "Can't create island set" );
         int startLoop = loopsindx[2 * loop_index];
         int endLoop   = loopsindx[2 * loop_index + 1];
  
         vector< EntityHandle > interiorCells;
         for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
         {
             EntityHandle cell = *cit;
             // see if it is in the interior of the loop
             CartVect center;
             MB_CHK_SET_ERR( mb->get_coords( &cell, 1, &( center[0] ) ), "Can't get cell center coords" );
             double lat = getLat( center ), lon = getLon( center );
             // if NA, use some boxes too, for lat/lon
             // if (NorthAmerica && (lat < 0.15 || lon < M_PI || lon > 2*M_PI - 0.5) )
             //  continue;
  
             if( interior_point( coords2d, startLoop, endLoop, lat, lon ) )
             {
                 interiorCells.push_back( cell );
                 MB_CHK_SET_ERR( mb->tag_set_data( tag1, &cell, 1, &loop_index ), "Can't get tag on cell" );
             }
         }
  
         MB_CHK_SET_ERR( mb->add_entities( islandSets[loop_index], &interiorCells[0], interiorCells.size() ),
                         "Can't add entities to set" );
     }
  
     std::stringstream islandFile;
  
     islandFile << "map.h5m";
  
     MB_CHK_SET_ERR( mb->write_file( islandFile.str().c_str() ), "Can't write island file" );
     return 0;
 }