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.

iphroPlacement.py

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()