Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
This is an implementation of the approach described in https://doi.org/10.1016/j.seta.2016.09.003 to determine good areas to place the IPHRO system. You'd need to download the GDEM data from NASA's website.
import numpy as np
import os
import math
import random
import overpy
from scipy.spatial import KDTree
def simulatedAnnealing(data, maxIter=1000000):
# NOTE: this function/approach wasn't effective enough at finding the optimal value
def pickRandomNeighbour(s):
i = random.randint(s[0] - 1, s[0] + 1)
j = random.randint(s[1] - 1, s[1] + 1)
while not (0 <= i < data.shape[0] and 0 <= j < data.shape[1] and (i, j) != s):
i = random.randint(s[0] - 1, s[0] + 1)
j = random.randint(s[1] - 1, s[1] + 1)
return (i, j)
def P(e, e_new, t):
if e_new > e:
return 1
return math.exp(-(e - e_new) / t)
s = (data.shape[0] // 2, data.shape[1] // 2)
for k in range(maxIter):
t = maxIter - k
s_new = pickRandomNeighbour(s)
if P(data[s[0], s[1]], data[s_new[0], s_new[1]], t) >= random.random():
s = s_new
return s
def approximateBestPlacement(data):
return np.unravel_index(np.argmax(data, axis=None), data.shape)
def getCoastalPoints():
api = overpy.Overpass()
res = api.query('relation["name:en"="Sinai Peninsula"];out geom;')
return np.array(
[
(float(geomValue.lat), float(geomValue.lon))
for member in res.relations[0].members
for geomValue in member.geometry
if float(geomValue.lat) < 29.4912 or float(geomValue.lat) > 31.0365
]
)
def getLatLong(i, j, baseLat, baseLong, sideLength):
"""
gets the lat/long points from a grid point
"""
unit = 1 / sideLength
return (baseLat + unit * (sideLength - i), baseLong + unit * j)
def getAIndex(
elevationData, coastalPointsKDTree: KDTree, baseLat, baseLong, sideLength
):
def getNearestMinDistance(point):
dd, _ = coastalPointsKDTree.query([point], k=1)
return dd[0]
def calcAIndex(elevation, distance):
if elevation < 500 or distance > 0.1:
return 0
return elevation / distance
return np.array(
[
[
calcAIndex(
elevationData[i, j],
getNearestMinDistance(
LatLongToXY(getLatLong(i, j, baseLat, baseLong, sideLength))
),
)
for j in range(elevationData.shape[1])
]
for i in range(elevationData.shape[0])
]
)
def LatLongToXY(
latlong,
):
return tuple(reversed(latlong))
def XYToLatLong(xy):
return tuple(reversed(xy))
def main():
coastalPoints = getCoastalPoints()
kdTree = KDTree([LatLongToXY(point) for point in coastalPoints])
for subdir, dirs, files in os.walk(".\GDEM_Data"):
for file in files:
baseLatitude = int(file[1:3]) if file[0] == "N" else -int(file[1:3])
baseLongitude = int(file[4:7]) if file[3] == "E" else -int(file[1:3])
name = os.path.join(subdir, file)
size = os.path.getsize(name)
dim = int(math.sqrt(size / 2)) # each int is two bytes
assert dim * dim * 2 == size
data = np.fromfile(name, np.dtype(">i2"), dim * dim).reshape((dim, dim))
aIndexes = getAIndex(
data, kdTree, baseLatitude, baseLongitude, data.shape[0]
)
np.savetxt(f"{file}_index.txt", aIndexes)
di, dj = approximateBestPlacement(aIndexes)
print(
f"{getLatLong(di, dj, baseLatitude, baseLongitude, data.shape[0])}, A-Index {aIndexes[di, dj]}, Elevation: {data[di, dj]}"
)
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment