Commit 1e33cdc4 authored by Fabian Kovac's avatar Fabian Kovac
Browse files

[i] general improvements

parent 885cdad8
......@@ -185,22 +185,22 @@ def get_distance(lon_a: np.array, lat_a: np.array, lon_b: np.array, lat_b: np.ar
def utm_to_lambert(lon: np.array, lat: np.array) -> tuple:
"""Convert UTM coordinates to a Lambert Conic Conformal Projection
"""Convert UTM coordinates to an Austrian Lambert Conic Conformal Projection (lcc)
EPSG code: 31287, see https://epsg.io/31287
Parameters:
lon (np.array): Vector containing longitudes
lat (np.array): Vector containing latitudes
lon (np.array): Vector containing utm longitudes
lat (np.array): Vector containing utm latitudes
Returns:
x, y (np.array): tuple of two vectors with coordinates
x, y (np.array): Vectors containing lcc coordinates
"""
# convert utm coordinates to radians
lon = np.radians(lon)
lat = np.radians(lat)
# define standard parallels according to EPSG:31287 - MGI/Austria Lambert
# --> see https://epsg.io/31287
lat1 = np.radians(49)
lat2 = np.radians(46)
......@@ -210,24 +210,28 @@ def utm_to_lambert(lon: np.array, lat: np.array) -> tuple:
lon0 = np.radians(13.33333333333333)
lat0 = np.radians(47.5)
# false easting/northing
# INCA grid: 701x401 km
# compensate for point of reference (middle of the grid)
# false easting: half of 701km = 350500m
# false northing: half of 401km = 200500m
# move point of reference (0, 0) to following x and y coordinates (x0, y0) to e.g. get strictly positive coordinates
# false easting: 3505000m
# false northing: 200500m
x0 = 350500
y0 = 200500
# volumetric mean radius of the earth in m
R = 6371000
# lambert conformal conic projection:
# --> see https://mathworld.wolfram.com/LambertConformalConicProjection.html
# lambert conformal conic projection as in https://pubs.usgs.gov/pp/1395/report.pdf
# reference:
# >> J. P. Snyder, Map projections--a working manual. Washington: U.S. G.P.O.,
# >> For sale by the Supt. of Docs., 1987., pp. 104-110
n = np.log(np.cos(lat1) * (1 / np.cos(lat2))) / np.log(np.tan(np.pi/4 + lat2/2) * (1 / np.tan(np.pi/4 + lat1/2)))
F = (np.cos(lat1) * np.tan(np.pi/4 + lat1/2)**n) / n
p = R * F * (1 / np.tan(np.pi/4 + lat/2))**n
p0 = R * F * (1 / np.tan(np.pi/4 + lat0/2))**n
# calculate x and y
# calculate x and y in lcc
x = p * np.sin(n * (lon - lon0)) + x0
y = p0 - p * np.cos(n * (lon - lon0)) + y0
......@@ -416,6 +420,10 @@ def prep() -> None:
_log('Calculated distances between sites using a WGS84 ellipsoid')
# convert CAPACITYINTERFACE to int
df_config['CAPACITYINTERFACE'] = df_config['CAPACITYINTERFACE'].map(lambda x: str(x)[:-5]).astype('int')
_log('Converted CAPACITYINTERFACE to int')
# convert FREQUENCY to float
df_config['FREQUENCY'] = df_config['FREQUENCY'].map(lambda x: str(x)[:-3]).astype('float')
_log('Converted FREQUENCY to float')
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment