diff --git a/.gitignore b/.gitignore index 55e94bb..465dc17 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ .env /dist *.cprof +*.pyc diff --git a/srtm/height_map_collection.py b/srtm/height_map_collection.py index 317c5bb..7f9c99c 100644 --- a/srtm/height_map_collection.py +++ b/srtm/height_map_collection.py @@ -94,23 +94,30 @@ def get_elevation_profile( end_latitude: float, end_longitude: float, apply_earth_curvature=True, + num_steps: int=10, ) -> List[ElevationProfilePoint]: """Get the elevation profile between the two points given""" values_per_degree = self.height_map_class.values_per_row - def to_int(lat_lng: float) -> int: - return round(lat_lng * values_per_degree) - - def to_float(lat_lng_int: int) -> float: - return lat_lng_int / values_per_degree - - points = points_on_line( - x1=to_int(start_latitude), - y1=to_int(start_longitude), - x2=to_int(end_latitude), - y2=to_int(end_longitude), - ) - converted_points = [(to_float(x), to_float(y)) for x, y in points] + # def to_int(lat_lng: float) -> int: + # return round(lat_lng * values_per_degree) + # + # def to_float(lat_lng_int: int) -> float: + # return lat_lng_int / values_per_degree + # + # points = points_on_line( + # x1=to_int(start_latitude), + # y1=to_int(start_longitude), + # x2=to_int(end_latitude), + # y2=to_int(end_longitude), + # ) + # converted_points = [(to_float(x), to_float(y)) for x, y in points] + + dx = float(end_longitude - start_longitude)/num_steps + dy = float(end_latitude - start_latitude)/num_steps + converted_points = [(start_latitude + _*dy, + start_longitude + _*dx) + for _ in range(num_steps + 1)] elevations = [] for latitude, longitude in converted_points: diff --git a/srtm/height_maps.py b/srtm/height_maps.py index 464f8b4..b416fdb 100644 --- a/srtm/height_maps.py +++ b/srtm/height_maps.py @@ -1,11 +1,45 @@ from pathlib import Path -from typing import Tuple, Callable +from typing import List, Tuple, Callable from zipfile import ZipFile +from math import floor from srtm.utilities import get_srtm3_file_path, get_srtm1_file_path from srtm.base_coordinates import RasterBaseCoordinates +# From: +def bilinear_interpolation(x: float, y: float, points: List[Tuple[float, float, float]]) -> float: + """ + Interpolate (x,y) from values associated with four points. + + The four points are a list of four triplets: (x, y, value). + The four points can be in any order. They should form a rectangle. + + >>> bilinear_interpolation(12, 5.5, + ... [(10, 4, 100), + ... (20, 4, 200), + ... (10, 6, 150), + ... (20, 6, 300)]) + 165.0 + + """ + # See formula at: http://en.wikipedia.org/wiki/Bilinear_interpolation + + points = sorted(points) # order points by x, then by y + (x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points + + if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2: + raise ValueError('points do not form a rectangle') + if not x1 <= x <= x2 or not y1 <= y <= y2: + raise ValueError('(x, y) not within the rectangle') + + return (q11 * (x2 - x) * (y2 - y) + + q21 * (x - x1) * (y2 - y) + + q12 * (x2 - x) * (y - y1) + + q22 * (x - x1) * (y - y1) + ) / ((x2 - x1) * (y2 - y1) + 0.0) + + class HeightMap: """Provides access to a single SRTM HGT file @@ -60,7 +94,7 @@ def validate(self): f"Expected {expected_bytes}, found {len(self.raster)}" ) - def get_altitude_for_pixel(self, x, y) -> int: + def get_altitude_for_pixel(self, x: int, y: int) -> int: """Get the height at the given pixel Will trigger loading of data @@ -68,7 +102,7 @@ def get_altitude_for_pixel(self, x, y) -> int: self.ensure_loaded() # Get the 1-indexed pixel number pixel_number = x + (y - 1) * self.values_per_row - # Convert it ro be 0-indexed + # Convert it to be 0-indexed pixel_number -= 1 byte_number = pixel_number * 2 @@ -76,12 +110,14 @@ def get_altitude_for_pixel(self, x, y) -> int: self.raster[byte_number : byte_number + 2], byteorder="big" ) - def get_altitude_for_latitude_and_longitude( - self, latitude: float, longitude: float - ) -> int: - """Get the height at the given lat/lng""" - x, y = self._latitude_and_longitude_to_coordinates(latitude, longitude) - return self.get_altitude_for_pixel(x, y) + # def get_altitude_for_latitude_and_longitude( + # self, latitude: float, longitude: float + # ) -> int: + # """Get the height at the given lat/lng""" + # x, y = self._latitude_and_longitude_to_coordinates(latitude, longitude) + # return self.get_altitude_for_pixel(x, y) + + _D_STEP = 0.0 def _latitude_and_longitude_to_coordinates( self, latitude: float, longitude: float @@ -89,8 +125,8 @@ def _latitude_and_longitude_to_coordinates( """Convert the given lat/long into x/y coordinates for this SRTM data""" origin_latitude = self.base_coordinates.latitude + 1 origin_longitude = self.base_coordinates.longitude - latitude_offset = origin_latitude - latitude - longitude_offset = longitude - origin_longitude + latitude_offset = origin_latitude - latitude + self._D_STEP*self.pixel_width + longitude_offset = longitude - origin_longitude + self._D_STEP*self.pixel_width if latitude_offset > 1 or latitude_offset < 0: raise ValueError( @@ -104,11 +140,41 @@ def _latitude_and_longitude_to_coordinates( ) # Add one because pixels are 1-indexed - x = round(longitude_offset / self.pixel_width) + 1 - y = round(latitude_offset / self.pixel_width) + 1 + x = floor(longitude_offset / self.pixel_width) + 1 + y = floor(latitude_offset / self.pixel_width) + 1 return x, y + def _coordinates_to_latitude_and_longitude( + self, x: int, y: int + ) -> Tuple[float, float]: + """Convert the given x/y coordinates into lat/long for this SRTM data""" + origin_latitude = self.base_coordinates.latitude + 1 + origin_longitude = self.base_coordinates.longitude + + longitude_offset = (x - 1)*self.pixel_width + latitude_offset = (y - 1)*self.pixel_width + + latitude = origin_latitude - latitude_offset + self._D_STEP*self.pixel_width + longitude = origin_longitude + longitude_offset - self._D_STEP*self.pixel_width + + return latitude, longitude + + def get_altitude_for_latitude_and_longitude( + self, latitude: float, longitude: float + ) -> int: + """Get the height at the given lat/lng""" + x, y = self._latitude_and_longitude_to_coordinates(latitude, longitude) + lat1, lon0 = self._coordinates_to_latitude_and_longitude(x, y) + lat0, lon1 = self._coordinates_to_latitude_and_longitude(x+1, y+1) + assert lat0 <= latitude <= lat1 + assert lon0 <= longitude <= lon1 + return bilinear_interpolation(latitude, longitude, + [(lat0, lon0, self.get_altitude_for_pixel(x, y+1)), + (lat1, lon0, self.get_altitude_for_pixel(x, y )), + (lat0, lon1, self.get_altitude_for_pixel(x+1, y+1)), + (lat1, lon1, self.get_altitude_for_pixel(x+1, y ))]) + class Srtm1HeightMap(HeightMap): """Provides access to a single SRTM HGT file