Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update high_places.py #1

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Update high_places.py #1

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
178 changes: 48 additions & 130 deletions high_places.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,138 +1,56 @@
#!/usr/bin/python3
##################################################################
# Tool to find high spots using 1deg x 1deg SRTM data
#
# The script will output a KML and a PNG file.
# The KML will display the terrain map with local maximum points
#
# Requirments = numpy, scipy, matplotlib, simplekml
#
# Source for SRTM data https://dwtkns.com/srtm30m/
#
#
###################################################################

import argparse
import numpy as np
from simplekml import Kml
from scipy.ndimage import maximum_filter
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource
import sys

# Sets the footprint of the local maximum box. Small numbers = more high points.
box_x = 300
box_y = 300

# Sets the color map of the terrain map
# See https://matplotlib.org/stable/tutorials/colors/colormaps.html
color_map = 'rainbow'

# I/O options
auto_name = True # Use this to automaticly name the output files based on input_filename

input_filename = 'N38W095.hgt' # ONLY USED IF auto_name=False SRTM File in .hgt format
## Note the output png and kml files will be auto named based on the input filename

print_list = False # Use this to toggle printing the list of high points to the command line

freedom_units = True # Use freedom (Imperical) units feet for elevation

###########################################################################################
###########################################################################################

# Each point in the 3600x3600 array = 1 deg / 3600 = 2.777777777777778e-4 deg
grid_step = 1/3600

# Filename as command line argument systax check
if len(sys.argv[1]) == 11:
if sys.argv[1][0] == 'N' or sys.argv[1][0] == 'S':
if sys.argv[1][3] == 'E' or sys.argv[1][3] == 'W':
input_filename = sys.argv[1]
else:
print('Filename failed syntax check')
quit()

# IMPORTANT the filename must be in N55E013.hgt repersenting North / South, 2 Digit Latitude in degrees, East / West, 3 Digit Longitude in degrees
print('Filename = %s converted to Lat = %d, Lon = %d' % (input_filename, int(input_filename[1:3]), int(input_filename[4:7])))

lat0 = int(input_filename[1:3])
lon0 = int(input_filename[4:7])

# North or South of equator, East or West of prime meridian
if input_filename[0] == 'N':
print('North Positive Lat')

if input_filename[0] == 'S':
print('South Negitive Lat')
lat0 = -abs(lat0)

if input_filename[3] == 'E':
print('East Positive Lon')

if input_filename[3] == 'W':
print('East Negitive Lon')
lon0 = -abs(lon0)

# Auto names output files if auto_name is True
image_filename = '%s.png' % (input_filename[0:7])
kml_filename = '%s.kml' % (input_filename[0:7])

# add 1 to lat0 due to hgt file origin being bottom left vs top right.
lat0_flipped = lat0 + 1

# Load elevation data
#elevations = np.loadtxt('N55E013.hgt', dtype='int16', skiprows=0)
elevations = np.fromfile(input_filename, np.dtype('>i2')).reshape((3601, 3601))
elevations = elevations.astype(float)

# Find local maxima using maximum_filter
footprint = np.ones((box_x, box_y))
local_max = (maximum_filter(elevations, footprint=footprint) == elevations) & (elevations > 0)

# Get indices of local maxima
max_indices = np.argwhere(local_max)

# Create a KML document
kml = Kml()

# Print number of hight points found
print('Found %d high points' % (len(max_indices)))

# generate terrain map you can change the color map
fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(elevations, cmap=color_map, vmin=elevations.min(), vmax=elevations.max())
plt.axis('off')
plt.savefig(image_filename, bbox_inches='tight', pad_inches=0, dpi=300)
plt.close()

print('Terrain map saved as %s' % (image_filename))

print('KML saved as %s' % (kml_filename))

# Generate the KML file
ground = kml.newgroundoverlay(name='Elevation Data')
ground.icon.href = image_filename
ground.gxlatlonquad.coords = [(lon0, lat0), (lon0 + 1, lat0), (lon0 + 1, lat0 + 1), (lon0, lat0 + 1)]

if print_list == True:
print('Lat, Lon, Elevation (m)')
if freedom_units == True:
print('Lat, Lon, Elevation (Ft)')

# Loop over indices and print corresponding elevation values
for index in max_indices:
row, col = index
elevation = elevations[row, col]
if freedom_units == True:
elevation = elevation * 3.28
lat = lat0_flipped - row * grid_step
lon = lon0 + col * grid_step
### Uncomment the next line for a verbose output
if print_list == True:
print('%1.4f, %1.4f, %d' % (lat, lon, elevation))
if freedom_units == False:
kml.newpoint(name=f'%dm' % (elevation), coords=[(lon, lat, elevation)])
if freedom_units == True:
kml.newpoint(name=f'%dft' % (elevation), coords=[(lon, lat, elevation)])

kml.save(kml_filename)
def parse_arguments():
parser = argparse.ArgumentParser(description='Tool to find high spots using 1deg x 1deg SRTM data.')
parser.add_argument('input_filename', help='SRTM File in .hgt format')
parser.add_argument('--box_x', type=int, default=300, help='Sets the footprint of the local maximum box. Small numbers = more high points.')
parser.add_argument('--box_y', type=int, default=300, help='Sets the footprint of the local maximum box. Small numbers = more high points.')
parser.add_argument('--color_map', default='rainbow', help='Sets the color map of the terrain map')
parser.add_argument('--print_list', action='store_true', help='Use this to toggle printing the list of high points to the command line')
parser.add_argument('--freedom_units', action='store_true', help='Use freedom (Imperical) units feet for elevation')
return parser.parse_args()

def load_elevation_data(filename):
elevations = np.fromfile(filename, np.dtype('>i2')).reshape((3601, 3601))
elevations = elevations.astype(float)
return elevations

def find_local_maxima(elevations, box_x, box_y):
footprint = np.ones((box_x, box_y))
local_max = (maximum_filter(elevations, footprint=footprint) == elevations) & (elevations > 0)
max_indices = np.argwhere(local_max)
return max_indices

def save_terrain_map(elevations, color_map, image_filename):
fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(elevations, cmap=color_map, vmin=elevations.min(), vmax=elevations.max())
plt.axis('off')
plt.savefig(image_filename, bbox_inches='tight', pad_inches=0, dpi=300)
plt.close()

def create_kml(elevations, max_indices, image_filename, kml_filename, lat0, lon0, freedom_units, print_list):
kml = Kml()
ground = kml.newgroundoverlay(name='Elevation Data')
ground.icon.href = image_filename
ground.gxlatlonquad.coords = [(lon0, lat0), (lon0 + 1, lat0), (lon0 + 1, lat0 + 1), (lon0, lat0 + 1)]

grid_step = 1 / 3600
lat0_flipped = lat0 + 1

for index in max_indices:
row, col = index
elevation = elevations[row, col]
if freedom_units:
elevation = elevation * 3.28
lat = lat0_flipped - row * grid_step
lon = lon0 + col * grid_step
if print_list:
print(f'{lat:.4f}, {lon:.4f}, {elevation:.0f}')
kml.newpoint(name=f'{elevation:.0f}{"ft" if freedom_units else "m"}', coords=[