Example scripts

On this page, some exemplary example scripts are presented, the scripts are shipped with the code.

Example 1 - Basic objects (GeoPoint, GeoVector3D)

Code

"""geonum example script 1

Introduction into the geonum base classes GeoPoint and GeoVector3D
"""

from geonum import GeoPoint
from SETTINGS import OPTPARSE
from numpy import testing as npt
from matplotlib.pyplot import show
### Create 2 GeoPoint objects 

# Summit region of Mt. Etna
p1 = GeoPoint(lat=37.751005, lon=14.993435, altitude=3264.0, name="Etna")

# Position of volcanological observatory of Mt. Etna
p2 = GeoPoint(37.765755,  15.016696, name="Observatory")

# Print info (string represenation of both points")
print(("Point1: %s, point 2: %s" %(p1, p2)))

# Get and print the connection vector of both points
connection_vector = p2 - p1

print(connection_vector)

# Import script options
(options, args) = OPTPARSE.parse_args()

# If applicable, do some tests. This is done only if TESTMODE is active: 
# testmode can be activated globally (see SETTINGS.py) or can also be 
# activated from the command line when executing the script using the 
# option --test 1
if int(options.test):
    from geonum import SRTM_AVAILABLE
    if not SRTM_AVAILABLE:
        raise ImportError('Skipping tests, since they require srtm.py to be installed')
    from os.path import basename
    
    actual = [connection_vector.azimuth,
              connection_vector.elevation, 
              connection_vector.magnitude]
    assert connection_vector.anchor is p1
    npt.assert_allclose(actual=actual,
                        desired=[51.378677249653983,
                                 -9.6064174085658465,
                                 2.6606074796318557],
                        rtol=1e-7)
    
    

    print(("All tests passed in script: %s" %basename(__file__))) 
try:
    if int(options.show) == 1:
        show()
except:
    print("Use option --show 1 if you want the plots to be displayed")



# Output:
# GeoVector3D Etna->Observatory
# Azimuth: 51.3786772497, Elevation: -9.60641740857, Magnitude: 2.66060747963 m


Example 2 - Introduction GeoSetup

Code

"""geonum example script 2

Introduction into the GeoSetup class
"""

from geonum import GeoPoint, GeoVector3D, GeoSetup, BASEMAP_AVAILABLE
from matplotlib.pyplot import show, close, rcParams
from os.path import join
from os import getcwd
from SETTINGS import OPTPARSE
from numpy import testing as npt
### Set save directory for figures
save_path = join(getcwd(), "scripts_out")

rcParams["font.size"] = 12
def create_geosetup():
    s = GeoSetup() 
    
    #Create two GeoPoints for source and instrument
    source = GeoPoint(37.751005,  14.993435, name="source")
    camera = GeoPoint(37.73122,  15.1129, name="cam")
    
    # Add the two GeoPoints to the GeoSetup
    s.add_geo_points(source, camera)
    s.set_borders_from_points()
    
    # Create plume vector anchored at source pointing south with horizontal
    # orientation (elevation angle 0)
    plume = GeoVector3D(azimuth=180, dist_hor=s.magnitude,
                        elevation=0, anchor=source, name="plume")
    # Create viewing direction vector anchored at instrument pointing west
    # at elevation angle of 8deg
    view_dir = GeoVector3D(azimuth=270, dist_hor=s.magnitude,
                           elevation=8, anchor=camera, name="cam_view")
            
    # Add the two GeoVectors to the GeoSetup class                        
    s.add_geo_vectors(plume, view_dir)
    return s
    
def plot_geosetup(geosetup):
    # Now plot 2D and 3D overview maps
    if not BASEMAP_AVAILABLE:
        raise ImportError('Basemap module is not available, cannot plot map')
    map2d = geosetup.plot_2d()
    map3d = geosetup.plot_3d()
    
    return map2d, map3d
    
if __name__ == "__main__":
    close("all")
    s = create_geosetup()
    try:
        map2d, map3d = plot_geosetup(s)
    except ImportError as e:
        print(repr(e))
    else:
        map2d.ax.figure.savefig(join(save_path, "ex2_out_1_map2D.png"))
        map3d.ax.figure.savefig(join(save_path, "ex2_out_2_map3D.png"))
        
    # Import script options
    (options, args) = OPTPARSE.parse_args()
    
    # If applicable, do some tests. This is done only if TESTMODE is active: 
    # testmode can be activated globally (see SETTINGS.py) or can also be 
    # activated from the command line when executing the script using the 
    # option --test 1
    if int(options.test):
        from geonum import SRTM_AVAILABLE
        if not SRTM_AVAILABLE:
            raise ImportError('Skipping tests, since they require srtm.py to be installed')
        from os.path import basename
        npt.assert_array_equal([4, 2, True],
                               [len(s.points),
                                len(s.vectors),
                                all([x in ['source', 
                                           'll', 
                                           'tr', 
                                           'cam'] for x in list(s.points)])],
                                all([x in ['plume', 
                                           'cam_view'] for x in list(s.vectors)]))
        
        actual = [s.points['source'].altitude,
                  s.points['cam'].altitude]
        npt.assert_allclose(actual=actual,
                            desired=[3264.0,
                                     803.0],
                            rtol=1e-7)
        print(("All tests passed in script: %s" %basename(__file__))) 
    try:
        if int(options.show) == 1:
            show()
    except:
        print("Use option --show 1 if you want the plots to be displayed")
    
    

Output figures

_images/ex2_out_1_map2D.png

Example 2, Fig. 1: 2D overview map of GeoSetup (Etna volcano)

_images/ex2_out_2_map3D.png

Example 2, Fig. 2: 3D overview map of GeoSetup (Etna volcano)

Example 3 - Mapping and plotting

Code

# -*- coding: utf-8 -*-
"""
geonum example script 3

Creating basemaps including topographic data and some further features 
(i.e. draw stuff into a map)
"""

import geonum
from numpy import mean
from matplotlib.pyplot import close, subplots, show
from os.path import join, exists
from os import getcwd
from SETTINGS import OPTPARSE
from numpy import testing as npt

### Set save directory for figures
save_path = join(getcwd(), "scripts_out")

close("all")
exceptions = []

def create_map_and_load_topodata():
    """Create and return map object around Guallatiri volcano"""
    # coordinate range of map
    lat0 = -18.48
    lon0 = -69.15
    lat1 = -18.39  
    lon1 = -69.05
    
    # trivial ...
    lat_center = mean([lat0, lat1])
    lon_center = mean([lon0, lon1])
    
    # Create a map object - the Map class is initialised as Basemap object
    # and is extended by some features
    m = geonum.mapping.Map(projection="lcc", llcrnrlat=lat0, 
                           llcrnrlon=lon0, urcrnrlat=lat1, urcrnrlon=lon1, 
                           lat_0=lat_center, lon_0=lon_center)
    
    # Compared to a normal Basemap object, a geonum Map can directly access 
    # and include SRTM topographic data
    m.load_topo_data("srtm")
    
    return m

def add_points_and_plot_map(basemap):
    """Add some points and plot in 2d and 3d"""
    basemap.draw_topo_3d()
    
    # create a geopoint for Guallatiri summit region 
    # (here, altitude is set manually)
    summit = geonum.GeoPoint(lat=-18.423672, lon=-69.090369,
                             altitude=6071.0, name="Guallatiri")
    
    # draw this point and a text into the map
    basemap.draw_geo_point_3d(summit)
    basemap.draw_text_3d(lon=-69.090369, lat=-18.423672, alt=6100.0,
                         text="Guallatiri", color="k")
    basemap.ax.set_title("Overview map Guallatiri volcano")
    
    # create two more objects (without specifying altitude -> is retrieved 
    # automatically from topo data)
    p1 = geonum.GeoPoint(-18.45, -69.12, name="Tourist (breathing hard)")
    p2 = geonum.GeoPoint(-18.40, -69.12, name="Llamas")
    
    
    # points can also be drawn directly into the map (here including the 
    # name which is placed 50 m above the topographic altitude of Tourist)
    p1.plot_3d(basemap, add_name=True, dz_text=50)
    p2.plot_3d(basemap, add_name=True, dz_text=50)
    
    # You can include a polygon connecting different points
    basemap.add_polygon_3d([summit, p2, p1], color="lime")
    
    # just some little adjustment of the 3D viewing direction
    basemap.ax.view_init(20, 240)
    
    
    # save the 3D figure
    print(("Save path: %s (exists: %s)" %(save_path, exists(save_path))))
    
    try:
        basemap.ax.figure.savefig(join(save_path, "ex3_out_1_map3D.png"))
    except Exception as e:
        exceptions.append("Failed to save first plot...%s" %repr(e))
        

    fig, ax = subplots(1,1)
    
    # update the axes object in the map
    basemap.ax = ax
    basemap.draw_topo() # draws topography into 2D map
    basemap.draw_topo_contour()
    
    basemap.draw_coordinates()
    basemap.draw_mapscale_auto()
    p1.plot_2d(basemap, add_name=True)
    p2.plot_2d(basemap, add_name=True)
    
    summit.plot_2d(basemap, add_name=True, dist_text=1.0, alpha=0.2,
                   angle_text=170, c="y")
    
    # You can include a polygon connecting different points (for whatever
    # reason)
    basemap.add_polygon_2d([summit, p2, p1], fc="lime", alpha=0.2)
    try:
        basemap.ax.figure.savefig(join(save_path, "ex3_out_2_map2D.png"))
    except Exception as e:
        exceptions.append("Failed to save second plot... %s" %repr(e))
    for e in exceptions:
        print(e)

if __name__ == "__main__":
    if not geonum.BASEMAP_AVAILABLE:
        raise ImportError('Cannot run script. Basemap is not installed')
    elif not geonum.SRTM_AVAILABLE:
        raise ImportError('Cannot run script. srtm.py is not installed')
    m = create_map_and_load_topodata()
    add_points_and_plot_map(m)
    
    # Import script options
    (options, args) = OPTPARSE.parse_args()
    
    # If applicable, do some tests. This is done only if TESTMODE is active: 
    # testmode can be activated globally (see SETTINGS.py) or can also be 
    # activated from the command line when executing the script using the 
    # option --test 1
    if int(options.test):
        from os.path import basename
        npt.assert_array_equal([m.topo_data.data.shape],
                               [(110, 122)])
        
        
        actual = [m.topo_data.data.mean(),
                  m.delta_lat,
                  m.delta_lon]
        
        actual.extend(m.get_map_center())
        npt.assert_allclose(actual=actual,
                            desired=[4982.759463487,
                                     0.089999999999,
                                     0.100000000000,
                                     -18.43500000000,
                                     -69.10000000000],
                            rtol=1e-7)
        print(("All tests passed in script: %s" %basename(__file__))) 
    try:
        if int(options.show) == 1:
            show()
    except:
        print("Use option --show 1 if you want the plots to be displayed")


Output figures

_images/ex3_out_1_map3D.png

Example 3, Fig. 1: Exemplary 3D overview map for Guallatiri volcano

_images/ex3_out_2_map2D.png

Example 3, Fig. 2: Exemplary 2D overview map for Guallatiri volcano

Example 4 - Elevation profiles and topographic analysis

Code

"""geonum example script 4

Elevation profiles and directional topographic analysis
"""

from matplotlib.pyplot import show, subplots
from os.path import join
from os import getcwd
from SETTINGS import OPTPARSE
from numpy import testing as npt

### Imports from other example scripts
from ex2_geosetup_intro import create_geosetup

### Set save directory for figures
save_path = join(getcwd(), "scripts_out")

def calc_profile():
    """Not very explanatory function name..."""
    
    
    # create Etna GeoSetup from example 2
    s = create_geosetup()
    
    # get camera
    cam = s.points["cam"]
    print(cam)
    
    # calculate elevation profile between camera and source (Etna summit)
    # the camera altitude is 15 m above topography
    elev_profile = cam.get_elevation_profile(geo_point=s.points["source"])
    
    
    return elev_profile
    
if __name__ == "__main__":
    profile = calc_profile()
    
    (dist, 
     dist_err, 
     intersect, 
     view_elevations, 
     ax)= profile.get_first_intersection(elev_angle=10.0,
                                         view_above_topo_m=15.0, 
                                         plot=True)
    
    ax.figure.savefig(join(save_path,
                           "ex4_out_1_elev_profile_intersect.png"))
    
    # Find the elvation angle corresponding to the horizon from the 
    # elevation profile
    elev, elev_secs, dist_secs = profile.find_horizon_elev(elev_start=10.0, 
                                                           elev_stop=20.0,
                                                           step_deg=0.1,
                                                           view_above_topo_m=15.0)
    
    fig, ax = subplots(1,1)
    ax.plot(elev_secs, dist_secs, "--x")
    ax.set_xlabel("Elevation angle [deg]")
    ax.set_ylabel("Retrieved distances [km]")
    ax.set_title("Elev horizon: %s" %elev)
    fig.savefig(join(save_path, "ex4_out_2_horizon_search.png"))
    
    # Import script options
    (options, args) = OPTPARSE.parse_args()
    
    # If applicable, do some tests. This is done only if TESTMODE is active: 
    # testmode can be activated globally (see SETTINGS.py) or can also be 
    # activated from the command line when executing the script using the 
    # option --test 1
    if int(options.test):
        from os.path import basename
        npt.assert_array_equal([], [])
        
        
        actual = [profile.start_point.longitude,
                  profile.start_point.latitude,
                  profile.start_point.altitude,
                  profile.alt_range,
                  dist, dist_err,
                  intersect.longitude,
                  intersect.latitude,
                  intersect.altitude,
                  elev]
        
        npt.assert_allclose(actual=actual,
                            desired=[15.1129, 37.73122, 803., 
                                     2482.5040909723493,
                                     7.667649e+00,   
                                     8.378121e-03,   
                                     1.502820e+01,   
                                     3.774504e+01,
                                     2.169563e+03,
                                     13.299999999999988],
                            rtol=1e-6)
        print(("All tests passed in script: %s" %basename(__file__))) 
    try:
        if int(options.show) == 1:
            show()
    except:
        print("Use option --show 1 if you want the plots to be displayed")
    

Output figures

_images/ex4_out_1_elev_profile_intersect.png

Example 4, Fig. 1: Exemplary elevation profile for GeoSetup from Ex.1 determined along azimuth angle connecting camera and summit including retrieval of intersection with topography for viewing elevation angle of 10 deg.

_images/ex4_out_2_horizon_search.png

Example 4, Fig. 2: Result from horizon search routine for profile shown in Fig. 1

Example 5 - Topography data access

Code

"""geonum example script 5

Low level access to topographic data (here SRTM) using the example of 
Oslo (Norway).

Note
----
Very short example, does not provide the deepest insights but gives an idea

"""
import geonum
from matplotlib.pyplot import show
from os.path import join
from os import getcwd
from SETTINGS import OPTPARSE
from numpy import testing as npt
### Set save directory for figures
save_path = join(getcwd(), "scripts_out")

def load_topo_oslo():
    dat = geonum.TopoDataAccess(mode='srtm')
    topodata = dat.get_data(59.8728, 10.375, 60.026, 10.9144)
    return topodata
    
    
if __name__ == "__main__":
    topo_data = load_topo_oslo()
    my_flat = geonum.GeoPoint(59.919386, 10.714970,
                              name="Somewhere in Frogner (Oslo)")
    
    if geonum.BASEMAP_AVAILABLE:
        basemap = topo_data.plot_2d()
        basemap.ax.set_title("SRTM topo data Oslo")
        basemap.draw_mapscale_auto()
        
        
        my_flat.plot_2d(basemap, add_name=True)
        basemap.ax.figure.savefig(join(save_path, "ex5_out_1_oslo_map.png"))
        
    # Import script options
    (options, args) = OPTPARSE.parse_args()
    
    # If applicable, do some tests. This is done only if TESTMODE is active: 
    # testmode can be activated globally (see SETTINGS.py) or can also be 
    # activated from the command line when executing the script using the 
    # option --test 1
    if int(options.test):
        from os.path import basename
        npt.assert_array_equal([], [])
        
        actual = [my_flat.altitude]
        npt.assert_allclose(actual=actual,
                            desired=[44.0],
                            rtol=1e-6)
        print(("All tests passed in script: %s" %basename(__file__))) 
    try:
        if int(options.show) == 1:
            show()
    except:
        print("Use option --show 1 if you want the plots to be displayed")
    
    

Output figures

_images/ex5_out_1_oslo_map.png

Example 5, Fig. 1: Overview map of Oslo area, Norway (SRTM dataset has some gaps here)