API Reference
Application programming interface of geonum.
Specifying locations on Earth
- class geonum.geopoint.GeoPoint(latitude=0, longitude=0, altitude=None, name=None, topo_access_mode=None, topo_path=None, topo_data=None, auto_topo_access=False)[source]
The Geopoint object represents a location in the atmosphere
This object is in 3D and includes altitude information.
- local_topo_path
directory that is used to search for local topography if topodata is requested. Irrelevant for default topo
- Type:
- Parameters:
latitude (float, optional) – latitude of point (decimal degrees), defaults to 0.
longitude (float, optional) – longitude of point (decimal degrees), defaults to 0.
altitude (float, optional) – elevation of point in m, defaults to 0.
name (str, optional) – name (ID) of this point, defaults to undefined
topo_access_mode (str, optional) – string specifying the current access mode for topographic data (in v1, choose between srtm or etopo1), defaults to None, in which case default of
TopoDataAccess
is used if topographic data is accessed via this point (e.g. if auto_topo_access is True). As of v1.4.X, default access mode is srtm, which uses SRTM database, which does not have full global coverage!topo_path (str, optional) – local path where local topography data is stored (at the moment only relevant for topo_mode etopo1), defaults to None.
topo_data (TopoData, optional) – existing topographic dataset that can be assigned to this point (can save time, e.g. for altitude access)
auto_topo_access (bool) – if True, try set
altitude`and :attr:`altitude_err
based on local topography (note that defa)
- offset(azimuth, dist_hor, dist_vert=0.0, ellipse=None, **kwargs) GeoPoint [source]
Returns new GeoPoint at offset position
- Parameters:
azimuth (float) – azimuth direction angle (in decimal degrees, e.g. 90 for E direction, -90 or 270 for west)
dist_hor (float) – horizontal offset to this point in km
dist_vert (float) – vertical offset to this point in m (positive if point is higher, negative, if it is lower)
ellipse (float) – geodetic system (ellipsoid), default is “WGS84”, i.e. World Geodetic System
**kwargs – additional keyword arguments passed to init of new
GeoPoint
object
- Returns:
new location at offset position
- Return type:
- set_topo_data(topo_data)[source]
Assign topographic data to this point
If input is valid, the topographic data set will be stored in class attribute topo_data.
- Parameters:
topo_data (TopoData) – topography dataset
- range_borders(*points, extend_fac=0.1)[source]
Get geographical borders (lower left, top right) of this point
Additional points can be included as non keyword arguments. The range borders are then determined under consideration of all specified points. If no additional points are specified, the range corners will be set corresponding to 1km diagonal distance to this point.
- get_topo_data(geo_point=None, azimuth=None, dist_hor=None, lon1=None, lat1=None)[source]
Retrieve topographic data grid spanned by this and another point
The second coordinate can be specified in several different ways dependent on which of the input arguments are used. First checks if the currently loaded topo data (
self.topo_data
) covers the range spanned by the 2 points. If this is the case, the loaded data will be used and nothing reloaded. If not, data will attempted to be loaded from the currentTopoDataAccess
object within the lon lat range covered by the 2 points (set usingrange_borders()
)Note
The following input combinations work (and are preferentially processed in the specified order if multiple input is given):
specify endpoint using geo_point
specify endpoint using azimuth and dist_hor
specify endpoint using lon1 and lat1
- Parameters:
geo_point (GeoPoint) – another geo_point,topo data will be retrieved between this point and destination point
azimuth (float) – azimuth angle of heading in decimal degrees (to be used with dist_hor)
dist_hor (float) – horizontal distance from this point in km (to be used with azimuth)
lon1 (float) – longitude of destination point, topo data will be retrieved between this point and destination point (to be used with lat1)
lat1 (float) – latitude of destination point, topo data will be retrieved between this point and destination point (to be used with lon1)
- Returns:
TopoData – topographic data object.
GeoPoint – the second coordinate which deterimes the domain wrt to this location.
- check_topo(lat1=None, lon1=None)[source]
Check if topography is available between this point and another
- get_elevation_profile(geo_point=None, azimuth=None, dist_hor=None, lon1=None, lat1=None, resolution=5.0, **mapping_opts)[source]
Estimates the elevation profile for a given viewing direction
The profile is retrieved up to a given distance from the coordinate of this point. For input possibilities see docs of
get_topo_data()
.Note
The following input combinations work (and are preferentially processed in the specified list order if multiple input is given):
specify endpoint using geo_point
specify endpoint using azimuth and dist_hor
specify endpoint using lon1 and lat1
- Parameters:
geo_point (GeoPoint) – end location of elevation profile (if this is provided, the following other input options to determine the end point will be ignored).
azimuth (float) – azimuth angle (direction) of elevation profile in decimal degrees (to be used with input arg dist_hor)
dist_hor (float) – horizontal distance from this point in km (to be used with azimuth)
lon1 (float) – longitude of destination point, topo data will be retrieved between this point and destination point (to be used with lat1)
lat1 (float) – latitude of destination point, topo data will be retrieved between this point and destination point (to be used with lon1)
resolution (float) – desired topo grid resolution in m. 1D interpolation of the elevation profile is performed if applicable.
- Returns:
the profile object
- Return type:
- set_topo_altitude()[source]
Set altitude using topographic terrain height at lat / lon position
The estimation is done by retrieving a 2x2 grid of topography data enclosing this point. The altitude value is estimated using the topo data tile closest to the coordinates of this point. The uncertainty is estimated conservatively using min/max difference of data in the 2x2 grid.
- Returns:
float – altitude
float – uncertainty in altitude value
- plot_2d(map, add_name=False, dist_text=0.5, angle_text=-45, **kwargs)[source]
Plot this point into existing 2D basemap
- Parameters:
- plot_3d(map, add_name=False, dz_text=0.0, **kwargs)[source]
Plot this point into existing 3D basemap
- mapbasemap
Basemap object (drawn in an Axes3D object)
- add_namebool
add the name of this GeoPoint in the map
- dz_textfloat
altitude offset of text (to point location in map)
- **kwargs
additional keyword arguments passed to matplotlib plot function.
Connecting the dots
- class geonum.geovector3d.GeoVector3D(dx=None, dy=None, dz=None, azimuth=None, dist_hor=None, elevation=None, anchor=None, name=None)[source]
A 3-dimensional geo vector object
3D vector representation for geodesic connections and arithmetics of locations (
geonum.GeoPoint
).Note
To create an instnace of
GeoVector3D
use one of the following to input combinations:dx, dy, dz
dz, azimuth, dist_hor
azimuth, dist_hor, elevation
Multiple input combinations will be processed in the preferred order as given in the list.
Note
Horitontal displacement components dx (longitude) and dy (latitude) are in units of km, while vertical displacement component dz (altitude) is in units of m.
- Parameters:
dx (float, optional) – longitudinal component in units of km
dy (float, optional) – latitudinal component in units of km
dz (float, optional) – altitude component in units of m
azimuth (float, optional) – azimuth orientation angle in units of decimal degrees, relative to North direction.
dist_hor (float, optional) – horizontal displacement length in input azimuth direction.
elevation (float, optional) – elevation angle in decimal degrees (0 corresponds to horizon, 90 to zenith)
anchor (GeoPoint, optional) – anchor point of this vector.
name (str, optional) – name of vector, defaults to “undefined”.
Example
>>> from geonum import GeoPoint, GeoVector3D >>> p = GeoPoint(latitude=10.0, longitude=15.0, name="random_point") >>> v = GeoVector3D(dx=15, dy=100, dz=300) #dx, dy in km, dz in m >>> new_point = p + v # GeoPoint object >>> print(new_point) GeoPoint undefined Lat: 10.904041412793307, Lon: 15.137202747069097, Alt: 300 m
- property norm: float
Norm of vector, wrapper for
magnitude()
- intersect_hor(other) GeoVector3D [source]
Find horizontal intersection of this vector with another vector
Note
Only works if anchor point (
anchor
) is set in both vectors.- Parameters:
other (GeoVector3D) – Other vector for which the intersection is to be determined.
- Raises:
AttributeError – if
anchor
is not set in this vector or input vector- Returns:
New vector pointing in direction of this vector but with correct horizontal magnitude to the intersection with the other vector.
- Return type:
Working with geographical domains
- class geonum.geosetup.GeoSetup(points=None, vectors=None, lat_ll=None, lon_ll=None, lat_tr=None, lon_tr=None, id=None, topo_access_mode=None, local_topo_path=None, cmap_vecs=None, init_borders=True)[source]
The GeoSetup class represents a collection of GeoPoints and vectors
- Parameters:
points (list) – list of
GeoPoint
objects to be included in this setupvectors (list) – list of
GeoVector3D
objects to be included in this setuplat_ll (float, optional) – lower left latitude of regime
lon_ll (float, optional) – lower left longitude of regime
lat_tr (float, optional) – top right latitude of regime
lon_tr (float, optional) – top right longitude of regime
id (str) – identification string of this setup
topo_access_mode (str) – topo data mode, default is SRTM (see
TopoDataAccess
for details)local_topo_path (str) – local path were topography data (e.g. ETOPO1 data) is stored
cmap_vecs (str) – String specifying a valid matplotlib colormap supposed to be used for drawing
GeoVector3D
objects into overview mapsinit_borders (bool) – Whether or not lower left (ll) and top right (tr) border points are supposed to be created.
- property magnitude
Returns dimension (in km) of area covered by this setup
- property topo_access
Topography data access class
- property cmap_vecs
Default colormap used for drawing vectors
- has_points()[source]
Determine whether any points are set in this GeoSetup
- Returns:
True if one or more GeoPoints are registered, else False
- Return type:
- has_topo_data()[source]
Check whether topographic data is assigned or not
- Returns:
True, if topo data is available, else not
- Return type:
- set_local_topo_path(p)[source]
Sets local path for Etopo1 data files can be found
Note
The default topo mode is “srtm” which provides online access, so it is not mandatory to provide topography data locally. However, the SRTM model has no global coverage, so there might be need to use another of the provided topo modes and provide the respective files locally.
- Parameters:
p (str) – new search path for topography data
- change_topo_mode(new_mode=None, local_path=None)[source]
Change the current mode for topography data access
- load_topo_data()[source]
Load topography data
Note
The loaded
TopoData
object will also be set in allGeoPoint
objects belonging to this setup
- add_geo_point(pt, assert_in_domain=True, overwrite_existing=False) None [source]
Add a GeoPoint to this collection
- Parameters:
pt (GeoPoint) – point to be added
assert_in_domain (bool) – if True, a check is performed whether the input point is within domain or not (does not apply if input point is “ll” or “tr”, that is, one of the points defining the domain itself). Defaults to True.
overwrite_existing (bool) – if True and a point with the same name already exists in this setup, then the existing point will be overwritten with input point, else an exception is raises. Defaults to False.
- Raises:
OutOfDomain – if paramerter assert_in_domain is True and input point is not within current domain.
ValueError – if parameter overwrite_existing is False and a point with the input n name already exists in this setup.
- Return type:
None
- add_geo_points(*pts, assert_in_domain=False) None [source]
Add multiple GeoPoints to the collection
- Parameters:
*pts – points to add
assert_in_domain (bool) – if True, check assert that each point is within domain
- Return type:
None
- add_geo_vector(vec, overwrite_existing=True) None [source]
Add
GeoVector3D
to this collection- Parameters:
vec (GeoVector3D) – vector to be added
overwrite_existing (bool) – If True and if vector with same name already exists, the former one
- Raises:
ValueError – if input is not GeoVector3D
- add_geo_vectors(*vecs) None [source]
Add multiple GeoVector3D objects to the collection
- Parameters:
*vecs – Instances of
GeoVector3D
to be added
- delete_geo_point(name) None [source]
Remove one of the geo points from the collection
- Parameters:
name (str) – name of the point
- Raises:
ValueError – if no point with with input name exists in this setup
- Return type:
None
- delete_geo_vector(name)[source]
Remove one of the vectors from the collection
- Parameters:
name (str) – name of the vector
- Raises:
ValueError – if no vector with with input name exists in this setup
- Return type:
None
- new_geo_point(*args, **kwargs) None [source]
Create new geo_point and add to collection
Create GeoPoint using input args and kwargs and then parse that GeoPoint to
add_geo_point()
.- Parameters:
*args – input args passed to init function of GeoPoint
**kwargs – input keyword args passed to init function of GeoPoint
- Return type:
None
- set_borders_from_points(extend_km=1, to_square=True) None [source]
Set lower left (ll) and top right (tr) corners of domain
The domain is inferred from all points associated with this setup.
- Parameters:
- Raises:
AttributeError – if no points are available in the setup.
- points_close(p, radius=None)[source]
Find all GeoPoints within a certain radius around input point
The search is performed against all points defined in this GeoSetup.
- Parameters:
- Returns:
names of all points that are in vicinity around input point
- Return type:
- static create_test_setup() GeoSetup [source]
Initiate example test data set
- Returns:
example setup
- Return type:
- plot_2d(draw_all_points=True, draw_all_vectors=True, draw_topo=True, draw_coastline=True, draw_mapscale=True, draw_legend=True, *args, **kwargs)[source]
To be updated
- plot_3d(draw_all_points=True, draw_all_vectors=True, cmap_topo='Oranges', contour_color='#708090', contour_lw=0.2, contour_antialiased=True, *args, **kwargs)[source]
Create a 3D overview map of the current setup
- Parameters:
draw_all_points (bool) – if True, all current GeoPoint objects are plotted, defaults to True
draw_all_vectors (bool) – if True, all current GeoVector3D objects are plotted, defaults to True
cmap (str) – string ID of the colormap used to plot the topographic data, defaults to “Oranges”
contour_color (str) – string specifying color of contour lines colors of contour lines (default: “#708090”)
contour_lw – width of drawn contour lines, defaults to 0.5, use 0 if you do not want contour lines inserted
contour_antialiased (bool) – apply antialiasing to surface plot of topography, defaults to False
*args – additional non-keyword parameters (passed to basemap)
**kwargs –
additional keyword parameters (passed to basemap)
- Returns:
plotted 3D basemap
- Return type:
Map
Computing elevation profiles
- class geonum.elevationprofile.ElevationProfile(observer, endpoint, topo_data=None, calc_on_init=True, **kwargs)[source]
Class for calculating elevation profiles
This class can be used to compute elevation profiles for arbitrary locations on earth. This is done by using topographic data which is provided (instance of
geonum.TopoData
). Profiles are computed with respect to input location (observer
) and the azimuthal direction of the profile is specified via the provided endpoint location (endpoint
).The profile is calculated between the two
GeoPoint
objects using a provided topographic data grid which must cover the range spanned by both points.Note
For class attributes see properties below.
- Parameters:
observer (GeoPoint) – start point of profile
endpoint (GeoPoint) – end point of profile
topo_data (TopoData) – topographic dataset
calc_on_init (bool, optional) – if True, then the profile is calculated on class inialisation. Default is True.
**kwargs – additional keyword args passed to
det_profile()
(only relevant if calc_on_init is True)
Example
>>> # define start and end location of profile >>> startpoint = GeoPoint(37.73122, 15.1129, name="Scientist") >>> endpoint = GeoPoint(37.751005, 14.993435, name="Etna Summit") >>> # Create profile retrieval engine >>> ep = ElevationProfile(observer=startpoint, endpoint=endpoint) >>> # Retrieve profile (uses default topographic data SRTM) >>> altitudes = ep.det_profile(interpolate=False) >>> # Print retrieved altitude levels along profile >>> print(altitudes)
- property topo_data: TopoData
Topographic data used to extract elevation profile
Note
Attempts automatic retrieval of topographic data if the data is not set.
- property line: LineOnGrid
- property profile: ndarray
Retrieved altitude levels in m along
line
- Raises:
AttributeError – if profile has not been calculated yet
- property dists: ndarray
Distances from observer in km along
line
- Raises:
AttributeError – if profile has not been calculated yet
- property resolution: float
Horizontal profile resolution (averaged from distance array)
Note
Only works if profile was already determined
- property slope: ndarray
Slope of profile at each point
The slope is calculated for each profile point by dividing the associated altitude step by the corresponding horitontal step.
- slope_angles(decimal_degrees=True)[source]
Get slope angles of profile (in each sample point)
- Parameters:
decimal_degrees (bool) – if True, angles are converted from radians to degrees
- Returns:
slope angles at each profile point
- Return type:
np.ndarray
- det_profile(interpolate=True, resolution=5.0, itp_type='linear', **mapping_opts)[source]
Determine the elevation profile
Searches the closest tiles in the topo data grid for both observer and endpoint and based on these two points the elevation profile is extracted using a
LineOnGrid
object (which extracts altitudes from the topo data along the connection vector of the 2 points using 2D spline intepolation)- Parameters:
interpolate (bool) – if True, the profile is interpolated to a certain horizontal resolution
resolution (float) – desired grid resolution in m for interpolation. Interpolation is performed if
interpolate
is True and if the actual resolution of the topo data is smaller than input, else, nothing is doneitp_type (str) – interpolation type (e.g. “linear”, “cubic”)
**mapping_opts – additional keyword arguments that are passed to
LineOnGrid.get_line_profile()
- Returns:
the array containing the retrieved altitude levels along the retrieval direction (from the observer)
- Return type:
np.ndarray
- get_altitudes_view_dir(elev_angle, view_above_topo_m=0)[source]
Get vector containing altitudes for a viewing direction
The viewing direction is specified by the azimuth angle of the connection vector between observer and endpoint, the elevation angle needs to be specified via the input parameters, and the start point (first elevation value) corresponds to the altitude at the observer position plus an offset in m which can be specified.
- get_first_intersection(elev_angle, view_above_topo_m=1.5, min_dist=None, local_tolerance=3, max_diff=None, plot=False)[source]
Finds first intersection of a viewing direction with topography
Start point of the viewing vector is the observer position (or more accurately, the center position of the closest topography tile in the topography data set). The relative altitude of the start point with respect to the topography altitude at the observer position can be controlled on input (arg
view_above_topo_m
) as well as the elevation angle of the viewing direction. The azimuth angle corresponds to the profile azimuth (i.e. azimuth betweenobserver
andendpoint
of this profile).The signal analysed for finding the intersection is a vector containing relative elevation positions of the viewing direction with respect to the profile altitude:
diff_signal = self.get_altitudes_view_dir - self.profile
The first intersection of the viewing direction with the topography is identified by the first zero crossing of this
diff_signal
.- Parameters:
elev_angle (float) – elevation angle of viewing direction (in decimal degrees)
view_above_topo_m (float) – altitude offset in m relative to altitude of start point (
observer
) in m. Defaults to 1.5 m.min_dist (float) – minimum distance (in km) of first considered intersection with topography from observer. If None, use 1% of distance between
observer
andendpoint
.local_tolerance (int) – tolerance factor to estimate distance uncertainty based on topo grid resolution. Defaults to 3.
max_diff (float, optional) – maximum allowed altitude difference in m between vector specifying viewing direction and the topography. Default to None, in which case the value assigned to
resolution
x 1000.plot (bool) – If true, the profile is plotted, including the intersection result.
- Returns:
float – retrieved horizontal distance of intersection point.
float – error of retrieved horizontal distance of intersection point.
GeoPoint – location of intersect
np.ndarray – elevations along viewing direction vector arising from input elevation and altitude offset as well as the azimuth of the profile
Axes or None – matplotlib axes instance if input arg plot is True, else None.
- find_horizon_elev(elev_start=0.0, elev_stop=60.0, step_deg=0.1, **kwargs)[source]
Find first elevation angle which does not intersect with topo
Scans towards zenith elevation within input elevation range in elevation steps provided via input arg step_deg, starting at elev_start. For each elevation, an attempt is made to find an intersection of the resulting viewing direction vector with the local topography using
get_first_intersection()
. The horizon elevation is identified as the first elevation angle for which no intersection can be found.- Parameters:
elev_start (float) – start search elevation angle
elev_stop (float) – stop search elevation angle
step_deg (float) – angle step for search (coarser is faster)
**kwargs – additional keyword agruments passed to
get_first_intersection()
- Raises:
IntersectNotFound – if no intersection with topography can be found for input elevation range.
- Returns:
float – detected elevation angle of horizon
list – list of elevation angles for which intersections with topography could be identified.
list – corresponding horizontal distances of
observer
to the terrain features.
- plot(ax=None, facecolor=None, alpha=0.2)[source]
Plot the elevation profile :param ax: matplotlib axes object to be used for plot :type ax: Axes, optional :param facecolor: color of profile, defaults to “#994d00” (light beige) :type facecolor: str, optional :param alpha: opaqueness of profile, defaults to 0.20. :type alpha: float, optional
- Returns:
matplotlib axes instance used for plot
- Return type:
Axes
Further tools
Line on grid
- class geonum.lineongrid.LineOnGrid(x0, y0, x1, y1, name=None)[source]
A class representing a line on a discrete grid
- Parameters:
Example
>>> import numpy as np >>> from geonum import LineOnGrid >>> data = np.ones((10,15)) * np.arange(15) >>> line = LineOnGrid(2,2,5,5) >>> profile = line.get_line_profile(data) >>> print(profile) [2. 2.74392912 3.50304939 4.24920976 5. ]
- property profile_coords: numpy.ndarray
(y,x) profile coordinates
- property length: int
Determine the length in units of grid indices
Note
The length is rounded to integer precision.
- get_line_profile(array, **kwargs)[source]
Retrieve line profile of data on a 2D array
Uses method
scipy.ndimage.map_coordinates()
Note
The spline interpolation will fail if there are NaNs in the input data. In this case, try using order=1 as additional input argument or use prefilter=False. For details see here
- Parameters:
array (numpy.ndarray) – 2D data array
**kwargs – additional keyword args that are passed to
scipy.ndimage.map_coordinates()
- Returns:
1D array containing the values along the line profile coordinates
- Return type:
- plot_line_on_grid(array, ax=None, **kwargs)[source]
Plot this line on an input array using imshow
- Parameters:
array (numpy.ndarray) – the data array to which this line is mapped
ax (matplotlib.axes.Axes, optional) – axes object
**kwargs – additional keyword arguments for imshow
- Returns:
axes instance used for plotting
- Return type:
- plot_line_profile(array, ax=None)[source]
Plot profile of line in input data array
- Parameters:
array (numpy.ndarray) – the data array to which this line is mapped
ax (matplotlib.axes.Axes, optional) – axes object
- Returns:
axes instance used for plotting
- Return type:
- plot(array)[source]
High level plot function, calls:
and puts them next to each other in a subplot
- Parameters:
array (numpy.ndarray) – 2D data array to which this line is mapped
- Returns:
figure instance containing the plots.
- Return type:
Topographic data
Access to topographic data
Access and handling of topographic data
- class geonum.topodataaccess.TopoDataAccess(mode=None, local_path=None)[source]
Factory class for accessing topographic data
This is a high-level factory class which handles the access of topography data. Registered topographic datasets can be found in
REGISTERED
Default access mode is SRTM (topographic dataset from NASA Shuttle Radar Topography Mission).
Example
>>> acc = TopoDataAccess(mode='srtm') >>> topo_data = acc.get_data(lat0=5, lon0=30, lat1=15, lon1=40)
Note
For developers: the registered access classes (such as
SRTMAccess
) should be based on (or follow the API of) template base classTopoAccessBase
.- mode
current access mode (string specifying which topographic dataset is supposed to be used for access).
- Type:
- Parameters:
mode (str) – one of the supported data access types (cf. keys of
REGISTERED
)local_path (str) – local path to etopo data (only relevant for etopo1 mode)
- REGISTERED = {'etopo1': <class 'geonum.topoaccessbase.Etopo1Access'>, 'srtm': <class 'geonum.topoaccessbase.SRTMAccess'>}
supported access modes (topographic datasets)
- DEFAULT_MODE = 'srtm'
- property modes
List of supported topographic datasets
- get_data(lat0, lon0, lat1=None, lon1=None, mode=None, local_path=None, check_allnan=True, **access_opts)[source]
Retrieve data from topography file
- Parameters:
lat0 (float) – start longitude for data extraction
lon0 (float) – start latitude for data extraction
lat1 (float) – stop longitude for data extraction (default: None). If None only data around lon0, lat0 will be extracted.
lon1 (float) – stop latitude for data extraction (default: None). If None only data around lon0, lat0 will be extracted
mode (str, optional) – mode specifying the topographic dataset that is supposed to be used
local_path (str, optional) – local path where topography data is stored (can be dictionary or filepath, is passed to corresponding access class and handled as implemented there)
check_allnan (bool) – if True and all retrieved values are NaN, then an error is thrown
**access_opts – additional access options that may be specific for the mode specified (e.g. search_database in case of etopo1)
- Returns:
object containing the data
- Return type:
- Raises:
TopoAccessError – if access fails
Topographic data class
- class geonum.topodata.TopoData(lats, lons, data, data_id=None, repl_nan_minval=False)[source]
Data class for topography data
This object represents topographic data in a certain latitude and longitude range, specified by the corner coordinates of the range (
lon0, lat0, lon1, lat1
). It may be used for 2D and 3D plotting of the topographic data (cf.plot_2d()
,plot_3d()
) or for further analysis such as coordinate based altitude retrievals (get_altitude()
), change of the grid resolution (cf.increase_grid_resolution()
) or for customised analysis by directly using the numpy data array containing the topographic altitudes (data
), together with the corresponding dimension coordinates which can be accessed via attributeslatitude`and :attr:`longitude
.To create an instance of this class, you may use the class
TopoDataAccess
.Note
This object does not provide topography data access. For data access, please use
TopoDataAccess
which returnsTopoData
objects. This object is intended for storage and post processing of topographic data- Parameters:
lats (ndarray) – numpy array with latitude coordinates of the topographic dataset. Accessible via
latitude
.lons (ndarray) – numpy array with longitude coordinates of the topographic dataset. Accessible via
longitude
.data (ndarray) – 2D numpy array containing elevation values, first index denotes latitude dimension, 2nd index denotes longitude dimension.
data_id (str) – ID of this data set.
repl_nan_minval (bool) – if True, then coordinates containing NaN values are replaced with the minimum altitude in data.
- property latitude
Wrapper for
lats
- property longitude
Wrapper for
lons
- replace_nans(fillval=None)[source]
Replace NaNs in topographic data with a fill value
- Parameters:
fillval (float) – value that is assigned to all coordinates that include NaNs
- property lon0
Returns first longitude (i.e.
self.lons[0]
)
- property lat0
Returns first latitude (i.e.
self.lats[0]
)
- property lon1
Returns last longitude (i.e.
self.lons[-1]
)
- property lat1
Returns last latitude (i.e.
self.lats[-1]
)
- property max
Returns maximum altitude of topo data
- property min
Returns minimum altitude of topo data
- property shape
Shape of 2D numpy array that contains topography
- property alt_range
Returns covered altitude range
- property resolution
Returns tuple (lat, lon) of the current grid resolution in km
Note
The resolution is determined at the center of this grid
- increase_grid_resolution(res=0.2, polyorder=2)[source]
Gaussian pyramide based upscaling of topographic grid
This function checks the current topographic resolution in the center of the grid. Based on this, an upsacling factor is determined and the
cv2.pyrUp()
is used to increase the resolution using interpolation. Note, that this does not increase the actual resolution of the topographic data grid.Note
Requires that opencv library is installed.
- Parameters:
- Returns:
new object with desired grid resolution
- Return type:
- Raises:
ImportError – if opencv is not installed.
- property delta_lon
Longitude range of data (in decimal degrees)
- property delta_lat
Latitude range of data (in decimal degrees)
- property center_coordinates
Tuple (lat, lon) with center coordinates of data
- init_mesh()[source]
Init X,Y meshgrid from latitudes and longitudes for map plots
- Returns:
meshgrid
- Return type:
- plot_3d(ax=None)[source]
Creates 3D surface plot of data
- Parameters:
ax – instance of matplotlib Axes3D object
- Returns:
plotted map object.
- Return type:
geonum.mapping.Map
- includes_coordinate(lat=None, lon=None)[source]
Checks if input coordinates are covered by this dataset
- closest_index(lat, lon)[source]
Finds the closest index to input coordinate
- Parameters:
- Returns:
2-element tuple containing the closest index of lat and lon arrays to input index
- Return type:
- Raises:
ValueError – if input coordinate is not included in this dataset
Low-level access and template classes
- class geonum.topoaccessbase.TopoAccessBase[source]
Abstract base class for topgraphy file implementations
Defines minimum interface for derived access classes of different topographic datasets.
- abstract get_data(lat0, lon0, lat1=None, lon1=None)[source]
Declaration of data access method
It is obligatory to implement this method into derived classes.
- Parameters:
lat0 (float) – first latitude coordinate of topographic range (lower left coord)
lon0 (float) – first longitude coordinate of topographic range (lower left coord)
lat1 (int or float, optional) – second latitude coordinate of topographic range (upper right coord). If None only data around lon0, lat0 will be extracted.
lon1 (int or float, optional) – second longitude coordinate of topographic range (upper right coord). If None only data around lon0, lat0 will be extracted.
- Returns:
instance of TopoData class
- Return type:
- class geonum.topoaccessbase.Etopo1Access(local_path=None, file_name=None, check_access=False, search_database=True)[source]
A class representing netCDF4 data access of Etopo1 data
See here for instructions on the data access.
- loader
data loader (
netCDF4.Dataset
)
- Parameters:
local_path (str) – directory where Etopo1 data files are stored
file_name (str) – file name of etopo data file
check_access (bool) – if True, then access to topography data is checked on init and an error is raised if no dataset can be accessed
search_database (bool) – if True and topodata file
file_path
does not exist, then a valid topography file is searched in all paths that are specified in file ~/.geonum/LOCAL_TOPO_PATHS.
- Raises:
TopoAccessError – if input arg check_access is True and if no supported data file can be found
- topo_id = 'etopo1'
ID of dataset
- supported_topo_files = ['ETOPO1_Ice_g_gmt4.grd', 'ETOPO1_Bed_g_gmt4.grd']
filenames of supported topographic datasets in preferred order
- property local_path
Directory containing ETOPO1 gridded data files
- property file_name
File name of topographic dataset used
- property file_path
Return full file path of current topography file
- get_data(lat0, lon0, lat1=None, lon1=None)[source]
Retrieve data from topography file
- Parameters:
lat0 (float) – first latitude coordinate of topographic range (lower left coord)
lon0 (float) – first longitude coordinate of topographic range (lower left coord)
lat1 (int or float, optional) – second latitude coordinate of topographic range (upper right coord). If None only data around lon0, lat0 will be extracted.
lon1 (int or float, optional) – second longitude coordinate of topographic range (upper right coord). If None only data around lon0, lat0 will be extracted.
- Returns:
instance of TopoData class
- Return type:
- class geonum.topoaccessbase.SRTMAccess(check_access=False, **kwargs)[source]
Class for SRTM topographic data access
Uses library srtm.py for online access of data.
Note
srtm.py
downloads the topo data from this source and stores a copy of the unzipped data files in the current cache directory found in home.Whenever data access is requested, the
srtm.py
checks if the file already exists on the local machine and if not downloads it online. The online access is rather slow, so do not be surprised, if things take a while when accessing a specific location for the first time.- Deleting cached SRTM files:
use
geonum.topoaccess.delete_all_local_srtm_files()
- Parameters:
check_access (bool) – check if data can be accessed on class initialisation
**kwargs – additional keyword arguments that are passed through (irrelevant for this class but relevant for factory loader class
TopoDataAccess
, particularlyset_mode()
therein.
- get_data(lat0, lon0, lat1=None, lon1=None)[source]
Load SRTM topographic subset for input range
- Parameters:
lat0 (float) – first latitude coordinate of topographic range (lower left coord)
lon0 (float) – first longitude coordinate of topographic range (lower left coord)
lat1 (int or float, optional) – second latitude coordinate of topographic range (upper right coord). If None only data around lon0, lat0 will be extracted.
lon1 (int or float, optional) – second longitude coordinate of topographic range (upper right coord). If None only data around lon0, lat0 will be extracted.
- Returns:
instance of TopoData class
- Return type:
Atmospheric calculations
- geonum.atmosphere.p0 = 101325.0
Pressure at sea level, Unit [Pa]
- geonum.atmosphere.NA = 6.022140857e+23
Avogadro constant, Unit [mol-1]
- geonum.atmosphere.T0_STD = 288.15
Standard sea level temperature, Unit [K]
- geonum.atmosphere.M_AIR_AVG = 28.9645
Average molar mass of air, Unit [g mol-1]
- geonum.atmosphere.R_STD = 8.31432
Gas constant (U.S. standard atmosphere), Unit [Nm mol-1 K-1]
- geonum.atmosphere.L_STD_ATM = -0.0065
Atmospheric lapse rate (Standard atmosphere), Unit [K m-1]
- geonum.atmosphere.L_DRY_AIR = -0.0098
Atmospheric lapse rate (dry atmosphere), Unit [K m-1]
- geonum.atmosphere.g(lat=45.0, alt=0.0)[source]
Gravitational accelaration as function of latitude and altitude
From Bodhaine et al., 1999, On Rayleigh Optical Depth Calculations who have it from List et al. 1968.
- Parameters:
- Returns:
Value(s) of Earth gravitational accelaration at sea level for input latitude(s) and altitude(s). If both inputs are arrays, then a 2D array is returned with the first axis (index 0) corresponding to altitudes and second axis to latitude
- Return type:
float or np.ndarray
Examples
>>> import geonum >>> import numpy as np >>> lats = np.linspace(0, 90, 200) >>> alts = np.linspace(0, 1000, 100) >>> dat = geonum.atmosphere.g(lats, alts) >>> dat.shape (100, 200)
- geonum.atmosphere.g0(lat=45.0)[source]
Gravitational accelaration at sealevel as function of latitude
See also
g()
for details- Parameters:
- Returns:
Value(s) of Earth gravitational accelaration at sea level for input latitude(s) and altitude(s). If both inputs are arrays, then a 2D array is returned with the first axis (index 0) corresponding to altitudes and second axis to latitude
- Return type:
float or np.ndarray
- geonum.atmosphere.temperature(alt=0.0, ref_temp=288.15, ref_alt=0.0, lapse_rate=-0.0065)[source]
Get temperature for a certain altitude (or altitude array)
Formula:
where:
is a reference temperature
is the atmospheric lapse-rate
is the altitude in m
is a reference altitude
- Parameters:
- Returns:
Temperature(s) in K corresponding to altitudes
- Return type:
float or np.ndarray
- geonum.atmosphere.beta_exp(mol_mass=28.9645, lapse_rate=-0.0065, lat=45.0)[source]
Exponent for conversion of atmosperic temperatures and pressures
Based on barometric height formula (see e.g. wikipedia) for details.
Formula:
where:
- geonum.atmosphere.pressure(alt=0.0, ref_p=101325.0, ref_temp=288.15, ref_alt=0.0, lapse_rate=-0.0065, mol_mass=28.9645, lat=45.0)[source]
Get atmospheric pressure in units of Pascal
Formula:
where:
is computed using
beta_exp()
is a reference pressure
is the temperature (cf.
temperature()
)is a reference temperature
is the atmospheric lapse-rate
is the altitude in m
is a reference altitude
- Parameters:
alt (float or np.ndarray) – altitude(s) in m
ref_p (float, optional) – reference pressure (default is std atm sea level)
ref_temp (float, optional) – reference temperature in K corresponding to ref_p (default is std atm sea level)
ref_alt (float, optional) – altitude corresponding to ref_p
lapse_rate (float, optional) – atmospheric lapse rate (in K m-1)
mol_mass (float, optional) – molar mass of air (in g mol-1)
lat (float, optional) – latitude in decimal degrees for calculation of gravitational constant
- Returns:
pressure(s) corresponding to input altitude(s)
- Return type:
float or np.ndarray
- geonum.atmosphere.pressure_hPa(alt=0.0, *args, **kwargs)[source]
Calls
pressure()
using provided input and returns in unit of hPaDefaults to standard atmosphere. This can be changed using further input parameters, for details see
pressure()
.- Parameters:
alt (float or np.ndarray) – altitude(s) in m
*args – non-keyword args parsed to
pressure()
**kwargs – keyword args parsed to
pressure()
- Returns:
pressure(s) in units of hPa corresponding to input altitude(s)
- Return type:
float or np.ndarray
- geonum.atmosphere.pressure2altitude(p, ref_p=101325.0, ref_temp=288.15, ref_alt=0.0, lapse_rate=-0.0065, mol_mass=28.9645, lat=45.0)[source]
General formula to convert atm. pressure to altitude
Formula:
where:
is a reference altitude
is a reference temperature
is the atmospheric lapse-rate
is the pressure (cf.
pressure()
)is a reference pressure
is computed using
beta_exp()
- Parameters:
p (float) – pressure in Pa
ref_p (float, optional) – reference pressure (default is std atm sea level)
ref_temp (float, optional) – reference temperature in K corresponding to ref_p (default is std atm sea level)
ref_alt (float, optional) – altitude corresponding to ref_p
lapse_rate (float, optional) – atmospheric lapse rate (in K / m)
mol_mass (float, optional) – molar mass of air
lat (float, optional) – latitude for calculation of gravitational constant
- Returns:
altitude corresponding to pressure level in defined atmosphere
- Return type:
Examples
>>> import geonum.atmosphere as atm >>> p = atm.pressure(2000) # pressure at 2000 m standard atm >>> int(atm.pressure2altitude(p)) 2000
- geonum.atmosphere.air_density(alt=0.0, temp=None, ref_p=101325.0, ref_temp=288.15, ref_alt=0.0, lapse_rate=-0.0065, mol_mass=28.9645, lat=45.0)[source]
Get atmospheric density in units of
Formula:
where:
is the atm. pressure (cf.
pressure()
)is the molar mass of air (defaults to standard atm)
is the gas constant (
R_STD
)is the temperature (cf.
temperature()
)
- Parameters:
alt (float or np.ndarray) – altitude(s) in m
temp (float, optional) – temperature in K, if None,
temperature()
is used to compute temperature.ref_p (float, optional) – reference pressure (default is std atm sea level)
ref_temp (float, optional) – reference temperature in K corresponding to ref_p (default is std atm sea level)
ref_alt (float, optional) – altitude corresponding to ref_p
lapse_rate (float, optional) – atmospheric lapse rate (in K / m)
mol_mass (float, optional) – molar mass of air
lat (float, optional) – latitude for calculation of gravitational constant
- Returns:
density corresponding to input altitude (in g m-3)
- Return type:
float or np.ndarray
- geonum.atmosphere.air_number_density(alt=0.0, temp=None, ref_p=101325.0, ref_temp=288.15, ref_alt=0.0, lapse_rate=-0.0065, mol_mass=28.9645, lat=45.0)[source]
Get atmospheric number density in m-3
- Parameters:
alt (float or np.ndarray) – altitude(s) in m
temp (float, optional) – temperature in K, if None,
temperature()
is used to compute temperature.ref_p (float, optional) – reference pressure (default is std atm sea level)
ref_temp (float, optional) – reference temperature in K corresponding to ref_p (default is std atm sea level)
ref_alt (float, optional) – altitude corresponding to ref_p
lapse_rate (float, optional) – atmospheric lapse rate (in K m-1)
mol_mass (float, optional) – molar mass of air
lat (float, optional) – latitude for calculation of gravitational constant
- Returns:
number density of air in #particles m-3 corresponding to altitude
- Return type:
float or np.ndarray
- geonum.atmosphere.refr_idx_300ppm_co2(lbda_mu=0.3)[source]
Get refractive index of air using eq. of Peck and Reeder, 1972
Uses either of two parametrisations for wavelengths above or below 0.23 microns using parametrisations from Peck and Reeder, 1972, Dispersion of air.
- geonum.atmosphere.refr_idx(lbda_mu=0.3, co2_ppm=400.0)[source]
Calculate refr. index of atm for varying CO2 amount
Uses either of two parametrisations for wavelengths above or below 0.23 microns from Peck and Reeder, 1972, Dispersion of air to retrieve the refractive index for CO2 concentrations at 300ppm and the parametrisation for varying CO2 amounts from Bodhaine et al., 1999, On Rayleigh Optical Depth Calculations.
- geonum.atmosphere.F_AIR(lbda_mu=0.3, co2_ppm=400.0)[source]
Depolarisation factor of air for Rayleigh scattering cross section
Taken from Bodhaine et al., 1999, On Rayleigh Optical Depth Calculations
- geonum.atmosphere.sigma_rayleigh(lbda_mu=0.3, co2_ppm=400.0)[source]
Rayleigh scattering cross section
Helper methods
- geonum.helpers.order_of_magnitude(num)[source]
Get exponent precision of input number
E.g. 1 if input number is 14.2, 0 if input number is 3 or -1 if input number is 0.1.
Note
The order of magnitude here is calculated using floor(log10(num)), thus, num=9, would result in a value of 0, in contrary to other definitions ( e.g. https://en.wikipedia.org/wiki/Order_of_magnitude).
- geonum.helpers.all_topodata_search_dirs()[source]
Returns a list of all directories that are searched for topography data
- Returns:
list containing all search directories (absolute paths)
- Return type:
- geonum.helpers.check_and_add_topodir(local_dir)[source]
Check if input directory is registered and if not, register it
- Parameters:
local_dir (str) – directory that is supposed to be checked
- geonum.helpers.isnum(val)[source]
Checks if input is number (int or float) or and not nan
- Parameters:
val – input object to be checked
- Returns:
whether or not input object val is of numerical type
- Return type:
- geonum.helpers.rotate_xtick_labels(ax, deg=30, ha='right')[source]
Rotate xtick labels in matplotlib axes object
- geonum.helpers.haversine_formula(lon0, lat0, lon1, lat1, radius=None)[source]
Haversine formula to compute distances on a sphere
Approximate horizontal distance between 2 points assuming a spherical earth.
- Parameters:
lon0 (float) – longitude of first point in decimal degrees
lat0 (float) – latitude of first point in decimal degrees
lon1 (float) – longitude of second point in decimal degrees
lat1 (float) – latitude of second point in decimal degrees
radius (float) – average earth radius in km, defaults to 6371 km
- Returns:
distance of both points in km
- Return type:
- geonum.helpers.approximate_connection_vector(lon0, lat0, lon1, lat1, len_lat_km=111.2)[source]
Returns approximate connection vector between two points
Note
Careful: Only approximate, suited for small distances ( < 100km)
- Parameters:
lon0 (float) – longitude of first point in decimal degrees
lat0 (float) – latitude of first point in decimal degrees
lon1 (float) – longitude of second point in decimal degrees
lat1 (float) – latitude of second point in decimal degrees
len_lat_km (float) – approx. length of 1 degree latitude at the equator in km.
- Returns:
2-element array containing (dx, dy) components of connection vector
- Return type:
- geonum.helpers.shifted_color_map(vmin, vmax, cmap=None)[source]
Shift center of a diverging colormap to value 0
Note
This method was found here (last access: 17/01/2017). Thanks to Paul H who provided it.
Function to offset the “center” of a colormap. Useful for data with a negative min and positive max and if you want the middle of the colormap’s dynamic range to be at zero level
- Parameters:
- Returns:
shifted colormap
- Return type:
Utility functions
Plotting
Helpers for plotting of maps etc
Note
BETA stage, code might undergo revisions
- geonum.plot_helpers.init_figure_with_geoaxes(projection=None, xlim=(-180, 180), ylim=(-90, 90))[source]
Initiate a matplotlib figure containing 1 cartopy GeoAxes instance
Can be used to plot maps.
- Parameters:
- Return type:
cartopy.mpl.geoaxes.GeoAxes
- geonum.plot_helpers.set_map_ticks(ax, xticks=None, yticks=None, tick_format=None)[source]
Set or update ticks in instance of GeoAxes object (cartopy)
Note
The input GeoAxes must be setup with
cartopy.crs.PlateCarree
projection.- Parameters:
ax (cartopy.GeoAxes) – map axes instance (needs to be set up with PlateCarree projection).
xticks (iterable, optional) – ticks of x-axis (longitudes). If None, it is retrieved automatically based on xlim of input axes. Defaults to None.
yticks (iterable, optional) – ticks of y-axis (latitudes). If None (or if xticks is None), it is retrieved automatically based on ylim of input axes. Defaults to None.
tick_formatter (str, optional) – string formatter for axes labels in the plot.
- Returns:
modified axes instance
- Return type:
cartopy.GeoAxes
- geonum.plot_helpers.plot_topo_contourf(ax, topo, oceans_separate=True, levels=50, cmap=None)[source]
- geonum.plot_helpers.plot_geopoint_into_map(ax, p, annotate=True, annot_kwargs=None, **kwargs)[source]
DEPRECATED: mapping.py
Note
The former map plotting routines routines (geonum.mapping
)
are deprecated as they are based on the basemap library. They
will be revised and replaced in parts using cartopy. For details see:
https://github.com/jgliss/geonum/issues/4