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.

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

Properly transform coordinates for geotiff contours #49

Merged
merged 2 commits into from
Jan 10, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
12 changes: 10 additions & 2 deletions pyhgtmap/__init__.py
Original file line number Diff line number Diff line change
@@ -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
Original file line number Diff line number Diff line change
@@ -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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Loading
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
Loading
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
Loading
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
Original file line number Diff line number Diff line change
@@ -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]'"
)
Loading
Loading