Skip to content

Commit

Permalink
Merge pull request #49 from agrenott/dev
Browse files Browse the repository at this point in the history
Properly transform coordinates for geotiff contours
  • Loading branch information
agrenott committed Jan 10, 2024
2 parents 89a2a83 + bee94cf commit b26f448
Show file tree
Hide file tree
Showing 17 changed files with 122 additions and 94 deletions.
12 changes: 10 additions & 2 deletions pyhgtmap/__init__.py
@@ -1,5 +1,5 @@
from importlib.metadata import version
from typing import List, Tuple
from typing import List, NamedTuple, Tuple

__author__ = "Aurélien Grenotton (agrenott@gmail.com)"
__version__ = version("pyhgtmap")
Expand All @@ -9,4 +9,12 @@
# Some type aliases
Polygon = List[Tuple[float, float]]
PolygonsList = List[Polygon]
BoudingBox = Tuple[float, float, float, float]


class BBox(NamedTuple):
"""Simple bounding box."""

min_lon: float
min_lat: float
max_lon: float
max_lat: float
13 changes: 6 additions & 7 deletions pyhgtmap/hgt/__init__.py
@@ -1,9 +1,8 @@
from __future__ import annotations

from typing import TYPE_CHECKING, Callable, Iterable, Tuple
from typing import Callable, Iterable, Tuple

if TYPE_CHECKING:
from pyhgtmap import BoudingBox
from pyhgtmap import BBox

# Coordinates transformation function prototype
TransformFunType = Callable[
Expand All @@ -12,7 +11,7 @@
]


def makeBBoxString(bbox: BoudingBox) -> str:
def makeBBoxString(bbox: BBox) -> str:
return f"{{0:s}}lon{bbox[0]:.2f}_{bbox[2]:.2f}lat{bbox[1]:.2f}_{bbox[3]:.2f}"


Expand All @@ -22,9 +21,9 @@ def transformLonLats(
maxLon: float,
maxLat: float,
transform: TransformFunType | None,
) -> BoudingBox:
) -> BBox:
if transform is None:
return minLon, minLat, maxLon, maxLat
return BBox(minLon, minLat, maxLon, maxLat)
else:
(lon1, lat1), (lon2, lat2), (lon3, lat3), (lon4, lat4) = transform(
[(minLon, minLat), (maxLon, maxLat), (minLon, maxLat), (maxLon, maxLat)],
Expand All @@ -33,4 +32,4 @@ def transformLonLats(
maxLon = max([lon1, lon2, lon3, lon4])
minLat = min([lat1, lat2, lat3, lat4])
maxLat = max([lat1, lat2, lat3, lat4])
return minLon, minLat, maxLon, maxLat
return BBox(minLon, minLat, maxLon, maxLat)
2 changes: 2 additions & 0 deletions pyhgtmap/hgt/contour.py
Expand Up @@ -130,6 +130,8 @@ def trace(self, elevation: int) -> tuple[list[numpy.ndarray], int, int]:
numOfPaths, numOfNodes = 0, 0
resultPaths = []
for path in rawPaths:
if self.transform:
path = numpy.array(self.transform(path))
path = simplify_path(path, self.rdp_epsilon)
splitPaths, numOfNodesAdd, numOfPathsAdd = self.splitList(path)
resultPaths.extend(splitPaths)
Expand Down
42 changes: 25 additions & 17 deletions pyhgtmap/hgt/file.py
Expand Up @@ -12,12 +12,13 @@
from matplotlib.path import Path as PolygonPath
from scipy import ndimage

from pyhgtmap import BBox
from pyhgtmap.hgt import TransformFunType, transformLonLats

from .tile import HgtTile

if TYPE_CHECKING:
from pyhgtmap import BoudingBox, Polygon, PolygonsList
from pyhgtmap import Polygon, PolygonsList
from pyhgtmap.cli import Configuration

with suppress(ImportError):
Expand Down Expand Up @@ -90,7 +91,7 @@ def parse_hgt_filename(
filename: str,
corrx: float,
corry: float,
) -> BoudingBox:
) -> BBox:
"""tries to extract borders from filename and returns them as a tuple
of floats:
(<min longitude>, <min latitude>, <max longitude>, <max latitude>)
Expand Down Expand Up @@ -123,7 +124,7 @@ def parse_hgt_filename(
f"something wrong with longitude coding in filename {filename:s}",
)
maxLon = minLon + 1
return minLon + corrx, minLat + corry, maxLon + corrx, maxLat + corry
return BBox(minLon + corrx, minLat + corry, maxLon + corrx, maxLat + corry)


def get_transform(
Expand Down Expand Up @@ -168,7 +169,7 @@ def parse_geotiff_bbox(
corrx: float,
corry: float,
doTransform: bool,
) -> BoudingBox:
) -> BBox:
try:
from osgeo import gdal, osr

Expand Down Expand Up @@ -210,7 +211,7 @@ def parse_geotiff_bbox(
maxLat,
transform,
)
return minLon + corrx, minLat + corry, maxLon + corrx, maxLat + corry
return BBox(minLon + corrx, minLat + corry, maxLon + corrx, maxLat + corry)
else:
# we need to take care for corrx, corry values then, which are always expected
# to be EPSG:4326, so transform, add corrections, and transform back to
Expand Down Expand Up @@ -239,15 +240,15 @@ def parse_geotiff_bbox(
maxLat,
reverseTransform,
)
return minLon, minLat, maxLon, maxLat
return BBox(minLon, minLat, maxLon, maxLat)


def parse_file_for_bbox(
fullFilename: str,
corrx: float,
corry: float,
doTransform: bool,
) -> BoudingBox:
) -> BBox:
fileExt: str = os.path.splitext(fullFilename)[1].lower().replace(".", "")
if fileExt == "hgt":
return parse_hgt_filename(os.path.split(fullFilename)[1], corrx, corry)
Expand All @@ -260,15 +261,15 @@ def calc_hgt_area(
filenames: list[tuple[str, bool]],
corrx: float,
corry: float,
) -> BoudingBox:
) -> BBox:
bboxes = [
parse_file_for_bbox(f[0], corrx, corry, doTransform=True) for f in filenames
]
minLon = sorted([b[0] for b in bboxes])[0]
minLat = sorted([b[1] for b in bboxes])[0]
maxLon = sorted([b[2] for b in bboxes])[-1]
maxLat = sorted([b[3] for b in bboxes])[-1]
return minLon, minLat, maxLon, maxLat
return BBox(minLon, minLat, maxLon, maxLat)


BBOX_EXPAND_EPSILON = 0.1
Expand Down Expand Up @@ -544,7 +545,7 @@ def init_as_geotiff(
else:
self.polygons = None

def borders(self, corrx=0.0, corry=0.0) -> BoudingBox:
def borders(self, corrx=0.0, corry=0.0) -> BBox:
"""determines the bounding box of self.filename using parseHgtFilename()."""
return parse_file_for_bbox(self.fullFilename, corrx, corry, doTransform=False)

Expand All @@ -558,7 +559,7 @@ def make_tiles(self, opts: Configuration) -> list[HgtTile]:

def truncate_data(
area: str | None, inputData: numpy.ma.masked_array
) -> tuple[BoudingBox, numpy.ma.masked_array]:
) -> tuple[BBox, numpy.ma.masked_array]:
"""truncates a numpy array.
returns (<min lon>, <min lat>, <max lon>, <max lat>) and an array of the
truncated height data.
Expand Down Expand Up @@ -628,12 +629,14 @@ def truncate_data(
maxLatTruncIndex:minLatTruncIndex,
minLonTruncIndex:maxLonTruncIndex,
]
return (realMinLon, realMinLat, realMaxLon, realMaxLat), zData
return BBox(realMinLon, realMinLat, realMaxLon, realMaxLat), zData
else:
return (self.minLon, self.minLat, self.maxLon, self.maxLat), inputData
return BBox(
self.minLon, self.minLat, self.maxLon, self.maxLat
), inputData

def chop_data(
inputBbox: BoudingBox,
inputBbox: BBox,
inputData: numpy.ma.masked_array,
depth=0,
):
Expand Down Expand Up @@ -666,7 +669,12 @@ def too_many_nodes(data: numpy.ma.masked_array) -> bool:
return False
return estim_num_of_nodes(data) > maxNodes

def get_chops(unchoppedData: numpy.ma.masked_array, unchoppedBbox):
def get_chops(
unchoppedData: numpy.ma.masked_array, unchoppedBbox
) -> tuple[
tuple[BBox, numpy.ma.masked_array],
tuple[BBox, numpy.ma.masked_array],
]:
"""returns a data chop and the according bbox. This function is
recursively called until all tiles are estimated to be small enough.
Expand All @@ -684,13 +692,13 @@ def get_chops(unchoppedData: numpy.ma.masked_array, unchoppedBbox):
unchoppedNumOfRows = unchoppedData.shape[0]
chopLatIndex = int(unchoppedNumOfRows / 2.0)
chopLat = unchoppedBboxMaxLat - (chopLatIndex * self.latIncrement)
lowerChopBbox = (
lowerChopBbox = BBox(
unchoppedBboxMinLon,
unchoppedBboxMinLat,
unchoppedBboxMaxLon,
chopLat,
)
upperChopBbox = (
upperChopBbox = BBox(
unchoppedBboxMinLon,
chopLat,
unchoppedBboxMaxLon,
Expand Down
6 changes: 3 additions & 3 deletions pyhgtmap/hgt/processor.py
Expand Up @@ -5,7 +5,7 @@
from multiprocessing.sharedctypes import Synchronized
from typing import TYPE_CHECKING, Callable, cast

from pyhgtmap import BoudingBox
from pyhgtmap import BBox
from pyhgtmap.hgt.file import HgtFile
from pyhgtmap.output.factory import get_osm_output

Expand Down Expand Up @@ -68,7 +68,7 @@ def single_output(self) -> bool:
def get_osm_output(
self,
hgt_files_names: list[str],
bounding_box: BoudingBox,
bounding_box: BBox,
) -> Output:
"""Allocate or return already existing OSM output (for consecutive calls in single output mode)
Expand Down Expand Up @@ -267,7 +267,7 @@ def process_files(self, files: list[tuple[str, bool]]) -> None:
self.get_osm_output(
[file_tuple[0] for file_tuple in files],
cast(
BoudingBox,
BBox,
[float(b) for b in self.options.area.split(":")],
),
)
Expand Down
9 changes: 5 additions & 4 deletions pyhgtmap/hgt/tile.py
Expand Up @@ -7,11 +7,12 @@
import numpy
import numpy.typing

from pyhgtmap import BBox
from pyhgtmap.hgt import TransformFunType, makeBBoxString, transformLonLats
from pyhgtmap.hgt.contour import ContoursGenerator, build_contours

if TYPE_CHECKING:
from pyhgtmap import BoudingBox, PolygonsList
from pyhgtmap import PolygonsList

meters2Feet = 1.0 / 0.3048

Expand All @@ -32,7 +33,7 @@ class HgtTile:

def __init__(
self,
bbox: BoudingBox,
bbox: BBox,
data: numpy.ma.masked_array,
increments: tuple[float, float],
polygons: PolygonsList | None,
Expand Down Expand Up @@ -88,7 +89,7 @@ def getElevRange(self) -> tuple[int, int]:
maxEle = int(self.zData.max())
return minEle, maxEle

def bbox(self, doTransform=True) -> BoudingBox:
def bbox(self, doTransform=True) -> BBox:
"""returns the bounding box of the current tile."""
if doTransform:
return transformLonLats(
Expand All @@ -99,7 +100,7 @@ def bbox(self, doTransform=True) -> BoudingBox:
self.transform,
)
else:
return self.minLon, self.minLat, self.maxLon, self.maxLat
return BBox(self.minLon, self.minLat, self.maxLon, self.maxLat)

def contourLines(
self,
Expand Down
6 changes: 3 additions & 3 deletions pyhgtmap/output/factory.py
Expand Up @@ -11,12 +11,12 @@
from . import make_elev_classifier

if TYPE_CHECKING:
from pyhgtmap import BoudingBox
from pyhgtmap import BBox
from pyhgtmap.cli import Configuration


def make_osm_filename(
borders: BoudingBox,
borders: BBox,
opts: Configuration,
input_files_names: list[str],
) -> str:
Expand Down Expand Up @@ -81,7 +81,7 @@ def makeBoundsString(bbox: Any) -> str:
def get_osm_output(
opts: Configuration,
input_files_names: list[str],
bounds: BoudingBox,
bounds: BBox,
) -> Output:
"""Return the proper OSM Output generator."""
outputFilename = make_osm_filename(bounds, opts, input_files_names)
Expand Down
4 changes: 2 additions & 2 deletions pyhgtmap/output/o5mUtil.py
Expand Up @@ -9,7 +9,7 @@
from pyhgtmap.varint import int2str, join, sint2str, writableInt, writableString

if TYPE_CHECKING:
from pyhgtmap import BoudingBox
from pyhgtmap import BBox
from pyhgtmap.hgt.tile import TileContours

HUNDREDNANO = 10000000
Expand Down Expand Up @@ -42,7 +42,7 @@ def __init__(
filename,
osmVersion,
pyhgtmap_version,
bbox: BoudingBox,
bbox: BBox,
elevClassifier: Callable[[int], str],
writeTimestamp=False,
) -> None:
Expand Down
4 changes: 2 additions & 2 deletions pyhgtmap/output/pbfUtil.py
Expand Up @@ -15,7 +15,7 @@
import pyhgtmap.output

if TYPE_CHECKING:
from pyhgtmap import BoudingBox
from pyhgtmap import BBox
from pyhgtmap.hgt.tile import TileContours


Expand All @@ -37,7 +37,7 @@ def __init__(
filename,
osmVersion,
pyhgtmap_version,
bbox: BoudingBox,
bbox: BBox,
elevClassifier: Callable[[int], str],
):
super().__init__()
Expand Down
Binary file modified tests/data/toulon_ref.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified tests/data/toulon_ref_smoothed.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added tests/data/toulon_ref_transformed.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
24 changes: 24 additions & 0 deletions tests/hgt/__init__.py
@@ -0,0 +1,24 @@
import contextlib
import importlib.util
from typing import Generator


@contextlib.contextmanager
def handle_optional_geotiff_support() -> Generator[None, None, None]:
"""
Context manager handling the cases where optional GeoTiff support has an impact.
Cases should run fully if geotiff dependencies are available, else specific exception is
expected.
"""
try:
# Execute test case
yield
except ImportError as ex:
if importlib.util.find_spec("osgeo") is not None:
# GDAL module is available, do NOT ignore the exception
raise
# GDAL not available, ensure the proper errror message is raised
assert ( # noqa: PT017 # Test is more complex
ex.msg
== "GeoTiff optional support not enabled; please install with 'pip install pyhgtmap[geotiff]'"
)

0 comments on commit b26f448

Please sign in to comment.