Orography / geometry

This submodule contains a class that converts between coordinates (wind-convention polar, cartesian, long/lat), and that handles the orography (query elevation relative to the base of the chimney at any coordinates).

`ensi_g14_calculator.orography`

Classes

Functions

Modules

`buildings`

Classes

`BuildingFootprint` `dataclass`

Rectangular building footprint with a flat roof, for Orography(..., buildings=[...]).

The centre is given either by latitude / longitude (degrees, like the chimneys and exposed groups) or by x / y (metres east / north of the stack; only meaningful for a single source). roof_elevation is absolute, in the datum of the terrain: metres above sea level for an elevation grid, the datum of flat_elevation for a flat orography.

    footprints = [
        BuildingFootprint(latitude=46.9301, longitude=7.4205, length=60 * u.m, width=25 * u.m,
                          orientation_deg=30.0, roof_elevation=560 * u.m),
    ]
    orography = Orography("terrain.asc", latitude0=..., longitude0=..., buildings=footprints)

Functions

centre_xy(orography)

Centre in the local frame of orography (metres east / north of its reference point).

contains(x_m, y_m, centre_m)

Boolean mask of the points (arrays in m, local frame) inside the footprint centred on centre_m.

Functions

`roof_elevation_m(footprints, centres_m, x_m, y_m)`

Highest roof elevation (m) of the footprints containing each point, NaN outside all footprints.

Parameters:

Name	Description	Default
`footprints`	sequence of BuildingFootprint	required
`centres_m`	their centres in the local frame, sequence of (x, y) in m	required
`x_m`	array of x coordinates (m)	required
`y_m`	array of y coordinates (m), same shape as x_m	required

`orography`

This module is concerned with the general geometry around our emission point.

That is:

orography (here, elevation in function of coordinates)
coordinates (lat/long to cartesian, cartesian to polar, and reverse)

Obtaining orography data

I suggest using the Shuttle Radar Topography Mission data, which has a permissive license. Its resolution is 30 meters, enough for this use case. Due to its nature (radar), it captures buildings and trees (not part of terrain), which for this use case is actually wanted.

To download the data, navigate to:

https://portal.opentopography.org/raster?opentopoID=OTSRTM.082015.4326.1

select the region of interest, and download as a .asc file. ("Arc ASCII Grid".)

Short usage example:

from ensi_g14_calculator.orography import Orography

orography = Orography(flat=True)  # for a flat orography
# for an orography with actual elevation data, you need to download a .asc file (see above) covering the
# region that interests you. Usually, 10 km x 10 km is more than enough. Then, instantiate your orography like this:
orography = Orography(asc_path="/path/to/downloaded/file.asc",
                      latitude0=latitude_of_release_point,
                      longitude0=longitude_of_release_point)

# then usually just pass the orography instance to other classes, although you can also call its
# methods directly if needed.

Classes

`Orography`

Class handling orography and coordinates (gives an elevation relative to a reference point), starting from an ascii file downloaded from opentopography and a reference point (longitude_0, latitude_0), typically the location of the chimney.

loads the grid
projects lon/lat -> local tangent-plane distances (x east, y north) from a reference point
interpolates elevation anywhere, either absolute or relative to the reference elevation

Flat-earth (equirectangular) approximation: good for up to 50 km scale.

Conventions:

longitude in degrees, latitude in degrees
x increases to the east, y increases to the north
x, y inputs should be astropy Quantities (e.g. m, m)
elevations are returned as floats in meters (same as the input grid)

If we do not wish to include detailed orographic data, we can pass flat=True instead and the resulting object will provide a flat orography. As a user, you should not have to use this class directly (it will be used by the other classes as the geometric reference). But you will need to instantiate it with your desired orography, and pass it to other classes.

Usage example:

orography = Orography(flat=True)  # for a flat orography
# for an orography with actual elevation data, you need to download a .asc file (see above) covering the
# region that interests you. Usually, 10 km x 10 km is more than enough. Then, instantiate your orography like this:
orography = Orography(asc_path="/path/to/downloaded/file.asc",
                      latitude0=latitude_of_release_point,
                      longitude0=longitude_of_release_point)

# then usually just pass the orography instance to other classes, although you can also call its
# methods directly if needed.

Functions

__init__(asc_path=None, latitude0=None, longitude0=None, *, origin='upper', check_reference=True, flat=False, flat_elevation=0.0 * u.m, buildings=None)

Parameters:

Name	Type	Description	Default
`buildings`		optional sequence of BuildingFootprint (orography/buildings.py). Inside a footprint the elevation is the absolute roof elevation, the terrain is only interpolated outside. The reference (stack base) stays the bare terrain at the stack, also for a stack on a roof: the stack height H_a remains the height above the ground, as in G-14.	`None`

latlon_to_xy(lat_deg, lon_deg)

Convert lon/lat (deg) -> (x,y) as Quantities in m.

Parameters:

Name	Type	Description	Default
`lat_deg`		latitude in degrees	required
`lon_deg`		longitude in degrees	required

Returns:

Type	Description
	tuple (x, y), cartesian coordinates in m with astropy units.

xy_to_latlon(x, y)

Convert (x,y) Quantities -> lon/lat (deg).

Parameters:

Name	Type	Description	Default
`x`		float or array, coordinate west->east, with astropy units of length	required
`y`		float or array, coordinate south_north, with astropy units of length	required

Returns:

Type	Description
	(latitude, longitude) tuple, in degrees.

elevation_xy(x, y, *, absolute_elevation=False)

Elevation at local (x,y) from reference.

Parameters:

Name	Description	Default
`x`	float or array, coordinate west->east, with astropy units of length	required
`y`	float or array, coordinate south_north, with astropy units of length	required
`absolute_elevation`	bool, default False. If true, returns the absolute elevation. If false, returns the elevation relative to the base of the stack of the chimney.	`False`

Returns:

Type	Description
	float or ndarray (meters), with astropy.unit attached (m)

elevation_latlon(lat_deg, lon_deg, *, absolute_elevation=False)

Elevation at (lon,lat) in degrees.

Parameters:

Name	Description	Default
`lat_deg`	latitude in degrees	required
`lon_deg`	longitude in degrees	required
`absolute_elevation`	bool, default False. If true, returns the absolute elevation. If false, returns the elevation relative to the base of the stack of the chimney.	`False`

Returns:

Type	Description
	float or ndarray (meters), with astropy.unit attached (m)

elevation_polar(angle, r, *, absolute_elevation=False, angle_unit='deg')

We match wind convention: 0° = North, 90° = East, clockwise positive

angle can be float/array. r can be Quantity/array.

Parameters:

Name	Description	Default
`angle`	polar angle in degrees, with 0° being north and 90° being east (clockwise -- matches wind conventions)	required
`r`	distance to stack, with astropy units (length)	required
`absolute_elevation`	bool, default False. If true, returns the absolute elevation. If false, returns the elevation relative to the base of the stack of the chimney.	`False`
`angle_unit`	'deg' or 'rad'	`'deg'`

Returns:

Type	Description
	float or ndarray (meters), with astropy.unit attached (m)

polar_to_xy(angle, r, angle_unit='deg')

We match wind convention: 0° = North, 90° = East, clockwise positive

Parameters:

Name	Description	Default
`angle`	polar angle, with 0° being north and 90° being east (clockwise -- matches wind conventions)	required
`r`	distance to stack, with astropy units (length)	required
`angle_unit`	string default 'deg', if 'deg', assumes angle given in degrees. if 'rad', assumes radians.	`'deg'`

Returns:

Type	Description
	x, y tuple, cartesian coordinates (with units of length)

xy_to_polar(x, y, angle_unit='deg')

Inverse of polar_to_xy using wind convention. Returns angle where 0° = North, 90° = East (clockwise).

Parameters:

Name	Description	Default
`x`	float or array, coordinate west->east, with astropy units of length	required
`y`	float or array, coordinate south_north, with astropy units of length	required
`angle_unit`	string, 'deg' or 'rad', unit of the angle.	`'deg'`

Returns:

Type	Description
	tuple (angle, r), the clockwise polar coordinates with 0° being north.

latlon_to_polar(lat_deg, lon_deg, angle_unit='deg')

Convert (lat,lon) to polar coordinates using wind convention. 0° = North, 90° = East, clockwise positive.

Parameters:

Name	Description	Default
`lat_deg`	latitude in degrees (scalar or array-like)	required
`lon_deg`	longitude in degrees (scalar or array-like), same shape as lat_deg	required
`angle_unit`	'deg' or 'rad', unit of the returned angle	`'deg'`

Returns:

Type	Description
	tuple (angle, r), where angle is clockwise from north and r is distance from reference

polar_to_latlon(angle, r, *, angle_unit='deg')

Convert polar coordinates to (lat,lon) using wind convention. 0° = North, 90° = East, clockwise positive.

Parameters:

Name	Description	Default
`angle`	polar angle in degrees or radians (scalar or array-like)	required
`r`	distance from reference point, with astropy units of length	required
`angle_unit`	'deg' or 'rad', unit of the input angle	`'deg'`

Returns:

Type	Description
	tuple (lat_deg, lon_deg), both in degrees

contains_latlon(lat_deg, lon_deg)

Quick check if a lon/lat point is inside the raster bounding box (not checking NODATA).

profile_along_line_xy(x0, y0, x1, y1, n=200, *, absolute_elevation=False)

Simple elevation profile between two local points. Returns (s, z) where s is distance along line (Quantity in m), z in meters.

Functions

`read_esri_ascii_grid(path)`

Loads an ESRI ASCII grid (.asc). This is apparently a fairly standard way of storing elevation maps, and as such our orography class should have the tool to load them.

Parameters:

Name	Type	Description	Default
`path`		string or Path, path to the location of the `.asc` file containing the orography data.	required

Returns:

Type	Description
	data (2D), extent (xmin_lon, xmax_lon, ymin_lat, ymax_lat), header (dict)

`plot_elevation(data, extent, ax=None, cmap='terrain', plot_hillshade=True)`

Just a diagnostic plot, plotting elevation data.

Parameters:

Name	Description	Default
`data`	2d array containing the elevation to be plotted	required
`extent`	as returned by the read_esri_ascii_grid function above, for the axes.	required
`ax`	the matplotlib axis to plot on, can be None (then one is created)	`None`
`cmap`	color map used for plotting	`'terrain'`
`plot_hillshade`	whether to add shades due to hills for added contrast.	`True`

Returns:

Type	Description
	None

Orography / geometry

ensi_g14_calculator.orography

Classes

Functions

Modules

buildings

Classes

BuildingFootprint dataclass

Functions

roof_elevation_m(footprints, centres_m, x_m, y_m)

orography

Obtaining orography data

Classes

Orography

Functions

read_esri_ascii_grid(path)

plot_elevation(data, extent, ax=None, cmap='terrain', plot_hillshade=True)

`ensi_g14_calculator.orography`

`buildings`

`BuildingFootprint` `dataclass`

`roof_elevation_m(footprints, centres_m, x_m, y_m)`

`orography`

`Orography`

`read_esri_ascii_grid(path)`

`plot_elevation(data, extent, ax=None, cmap='terrain', plot_hillshade=True)`