Update download of elevation data + other fixes

requirements.txt

@@ -1,10 +1,8 @@
-certifi==2023.7.22
-chardet==4.0.0
-idna==2.10
-numpy==1.21.0
-Pillow==10.2.0
-pyproj==3.1.0
-requests==2.31.0
-Rtree==0.9.7
-Shapely==1.7.1
-urllib3==1.26.18
+numpy==2.2.6
+Pillow
+pyproj
+requests
+Rtree
+Shapely
+scikit-fmm
+rasterio

src/O4_DEM_Utils.py

@@ -7,13 +7,7 @@
 from math import sqrt
 import array
 import numpy
-
-try:
-    from osgeo import gdal
-    has_gdal = True
-    gdal.UseExceptions()
-except:
-    has_gdal = False
+import rasterio
 from PIL import Image
 import O4_UI_Utils as UI
 import O4_File_Names as FNAMES
@@ -501,31 +495,35 @@ def read_elevation_from_file(
         x1 = y1 = 1
         epsg = 4326
         nodata = -32768
-    elif has_gdal:
+    else:
         try:
-            ds = gdal.Open(file_name)
-            rs = ds.GetRasterBand(1)
+            dataset = rasterio.open(file_name)
             if not info_only:
-                alt_dem = rs.ReadAsArray().astype(numpy.float32)
-            (nxdem, nydem) = (ds.RasterXSize, ds.RasterYSize)
-            nodata = rs.GetNoDataValue()
+                alt_dem= dataset.read(1).astype(numpy.float32)
+            (nxdem, nydem)=(dataset.width, dataset.height)
+            nodata=dataset.nodata
             if nodata is None:
                 UI.vprint(
                     1,
                     "    WARNING: raster DEM does not advertise its no_data ",
                     "value, assuming -32768.",
                 )
                 nodata = -32768
-            else:  
-                # elevations being stored as float32, we push the nodata to that 
-                # framework too, and then replace no_data values by -32768 
+            else:
+                # elevations being stored as float32, we push the nodata to that
+                # framework too, and then replace no_data values by -32768
                 # anyway for uniformity
                 nodata = numpy.float32(nodata)
                 if not info_only:
                     alt_dem[alt_dem == nodata] = -32768
                 nodata = -32768
+
+            #Replace nan with nodata
+            alt_dem=numpy.nan_to_num(alt_dem,nan=nodata)
+
             try:
-                epsg = int(ds.GetProjection().split('"')[-2])
+                crs=dataset.crs
+                epsg=crs.to_epsg()
             except:
                 UI.vprint(
                     1,
@@ -536,8 +534,8 @@ def read_elevation_from_file(
             if epsg not in (
                 4326,
                 4269,
-            ):  
-            # let's be blind about 4269 which might be sufficiently close to 
+            ):
+            # let's be blind about 4269 which might be sufficiently close to
             # 4326 for our purposes
                 UI.lvprint(
                     1,
@@ -546,7 +544,8 @@ def read_elevation_from_file(
                     ". Only EPSG:4326 is supported, result is likely to ",
                     "be non sense.",
                 )
-            geo = ds.GetGeoTransform()
+            transform=dataset.transform
+            geo=(transform.c, transform.a, transform.b, transform.f, transform.d, transform.e)
             # We are assuming AREA_OR_POINT is area here
             x0 = geo[0] + 0.5 * geo[1] - lon
             y1 = geo[3] + 0.5 * geo[5] - lat
@@ -568,22 +567,6 @@ def read_elevation_from_file(
             x1 = y1 = 1
             epsg = 4326
             nodata = -32768
-    elif not has_gdal:
-        UI.lvprint(
-            1,
-            "   WARNING: unsupported raster (install Gdal):",
-            file_name,
-            "-> replaced with zero altitude.",
-        )
-        nxdem = nydem = base_if_error
-        if not info_only:
-            alt_dem = numpy.zeros(
-                (base_if_error, base_if_error), dtype=numpy.float32
-            )
-        x0 = y0 = 0
-        x1 = y1 = 1
-        epsg = 4326
-        nodata = -32768
     return (epsg, x0, y0, x1, y1, nodata, nxdem, nydem, alt_dem)
 
 
@@ -592,90 +575,64 @@ def read_elevation_from_file(
 ##############################################################################
 def ensure_elevation(source, lat, lon, verbose=True):
     if source == "View":
-        # Viewfinderpanorama grouping of files and resolutions is a 
+        # Viewfinderpanorama grouping of files and resolutions is a
         # bit complicated...
-        if (lat, lon) in (
-            (44, 5),
-            (45, 5),
-            (46, 5),
-            (43, 6),
-            (44, 6),
-            (45, 6),
-            (46, 6),
-            (47, 6),
-            (43, 7),
-            (44, 7),
-            (45, 7),
-            (46, 7),
-            (47, 7),
-            (45, 8),
-            (46, 8),
-            (47, 8),
-            (45, 9),
-            (46, 9),
-            (47, 9),
-            (45, 10),
-            (46, 10),
-            (47, 10),
-            (45, 11),
-            (46, 11),
-            (47, 11),
-            (45, 12),
-            (46, 12),
-            (47, 12),
-            (46, 13),
-            (47, 13),
-            (46, 14),
-            (47, 14),
-            (46, 15),
-            (47, 15),
-        ):
-            resol = 1
-            url = (
-                "http://viewfinderpanoramas.org/dem1/"
-                + os.path.basename(FNAMES.base_file_name(lat, lon)).lower()
-                + ".zip"
-            )
+
+        deferranti_nbr = 31 + lon // 6
+        if deferranti_nbr < 10:
+            deferranti_nbr = "0" + str(deferranti_nbr)
         else:
-            deferranti_nbr = 31 + lon // 6
-            if deferranti_nbr < 10:
-                deferranti_nbr = "0" + str(deferranti_nbr)
-            else:
-                deferranti_nbr = str(deferranti_nbr)
-            alphabet = list("ABCDEFGHIJKLMNOPQRSTUVWXYZ")
-            deferranti_letter = (
-                alphabet[lat // 4] if lat >= 0 else alphabet[(-1 - lat) // 4]
-            )
-            if lat < 0:
-                deferranti_letter = "S" + deferranti_letter
-            if deferranti_letter + deferranti_nbr in (
-                "O31",
-                "P31",
-                "N32",
-                "O32",
-                "P32",
-                "Q32",
-                "N33",
-                "O33",
-                "P33",
-                "Q33",
-                "R33",
-                "O34",
-                "P34",
-                "Q34",
-                "R34",
-                "O35",
-                "P35",
-                "Q35",
-                "R35",
-                "P36",
-                "Q36",
-                "R36",
+            deferranti_nbr = str(deferranti_nbr)
+        alphabet = list("ABCDEFGHIJKLMNOPQRSTUVWXYZ")
+        deferranti_letter = (
+            alphabet[lat // 4] if lat >= 0 else alphabet[(-1 - lat) // 4]
+        )
+        if lat < 0:
+            deferranti_letter = "S" + deferranti_letter
+        if deferranti_letter + deferranti_nbr in (
+               # "L31",
+               # "L32",
+               # "L33",
+               # "K32",
+               # "O31",
+               # "P31",
+               # "N32",
+               # "O32",
+               # "P32",
+               # "Q32",
+               # "N33",
+               # "O33",
+               # "P33",
+               # "Q33",
+               # "R33",
+               # "O34",
+               # "P34",
+               # "Q34",
+               # "R34",
+               # "O35",
+               # "P35",
+               # "Q35",
+               # "R35",
+               # "P36",
+               # "Q36",
+               # "R36",
+                 #Greenland
+                "U19","U20","U21","U22","U23","U24","U25","U26","U27","U28","U29",
+                "T18","T19","T20","T21","T22","T23","T24","T25","T26","T27","T28",
+                "S19","S20","S21","S22","S23","S24","S25","S26","S27","S28",
+                "R21","R22","R23","R24","R25","R26","R27",
+                "Q22","Q23","Q24","Q25",
+                "P22","P23","P24",
+                "O23",
             ):
                 resol = 1
-            else:
+        else:
                 resol = 3
-            url = (
+                # 3" DEM elevation data for Iceland
+                if deferranti_letter+str(deferranti_nbr) in ("Q27","Q28"):
+                   deferranti_letter = "ISL"
+                   deferranti_nbr = ""
+        url = (
                 "http://viewfinderpanoramas.org/dem"
                 + str(resol)
                 + "/"
@@ -719,7 +676,7 @@ def ensure_elevation(source, lat, lon, verbose=True):
                 if ("W" in fname) or ("w" in fname):
                     lon0 *= -1
                 out_filename = FNAMES.viewfinderpanorama(lat0, lon0)
-                # we don't wish to overwrite a 1" version by downloading 
+                # we don't wish to overwrite a 1" version by downloading
                 # the whole archive of a nearby 3" one
                 if (
                     not os.path.exists(out_filename)
@@ -741,7 +698,7 @@ def ensure_elevation(source, lat, lon, verbose=True):
             "    WARNING : This elevation source has no longer direct downloads !"
         )
         return 0
-        # TODO : is there a way to get it back (worth it ?) 
+        # TODO : is there a way to get it back (worth it ?)
         url = "https://cloud.sdsc.edu/v1/AUTH_opentopography/Raster/"
         if source == "SRTM":
             url += "SRTM_GL1/SRTM_GL1_srtm/"

src/O4_Mask_Utils.py

@@ -172,7 +172,7 @@ def build_mask(til_x, til_y):
         )[1024 : 4096 + 1024, 1024 : 4096 + 1024]
 
         if tile.masks_custom_extent:
-            blured_mask = numpy.maximum(blured_mask, custom_mask)
+            blured_mask = numpy.maximum(blured_mask, custom_array)
 
         if not (blured_mask.max() == 0 or blured_mask.min() == 255):
             mask_im = Image.fromarray(blured_mask)

src/O4_Mesh_Utils.py

@@ -689,6 +689,7 @@ def build_mesh(tile):
     fingers_crossed = subprocess.Popen(
         mesh_cmd, stdout=subprocess.PIPE, bufsize=0
     )
+
     while True:
         line = fingers_crossed.stdout.readline()
         if not line:
@@ -700,29 +701,43 @@ def build_mesh(tile):
                 pass
     time.sleep(0.3)
     fingers_crossed.poll()
+
     if fingers_crossed.returncode:
-        UI.vprint(
-            0,
-            "\nWARNING: Triangle4XP could not achieve the requested quality ",
-            "(min_angle), most probably due to an uncatched OSM error.\n",
-            "It will be tempted now with no angle constraint ",
-            "(i.e. min_angle=0).",
-        )
-        mesh_cmd[-5] = "{:.9g}".format(0)
-        fingers_crossed = subprocess.Popen(
-            mesh_cmd, stdout=subprocess.PIPE, bufsize=0
-        )
-        while True:
-            line = fingers_crossed.stdout.readline()
-            if not line:
+        min_angles_to_try = [5, 2, 0]
+        for mina in min_angles_to_try:
+            print(
+                "\n\nWARNING: Triangle4XP could not achieve the requested quality",
+                "(min_angle), most probably due to an uncatched OSM error\n",
+                f"        It will be retried with lower value (min_angle = {mina})\n\n",
+            )
+
+            Tri_option = (
+               "-pq" + "{:.9g}".format(mina) + do_refine +
+               "uYB" + tri_verbosity + output_poly + limit_tris
+            )
+
+            mesh_cmd[1] = Tri_option
+
+            UI.vprint(2, "   Mesh command:", " ".join(mesh_cmd))
+            fingers_crossed = subprocess.Popen(mesh_cmd, stdout=subprocess.PIPE, bufsize=0)
+
+            while True:
+                line = fingers_crossed.stdout.readline()
+                if not line:
+                    break
+                else:
+                    try:
+                        print(line.decode("utf-8")[:-1])
+                    except:
+                        pass
+            time.sleep(0.3)
+            fingers_crossed.poll()
+
+            if fingers_crossed.returncode == 0:
+                # If return code is 0, process succeeded
+                print(f"\nSuccess with min_angle = {mina}\n")
                 break
-            else:
-                try:
-                    print(line.decode("utf-8")[:-1])
-                except:
-                    pass
-        time.sleep(0.3)
-        fingers_crossed.poll()
+
         if fingers_crossed.returncode:
             UI.exit_message_and_bottom_line(
                 "\nERROR: Triangle4XP really couldn't make it !\n\n",

src/O4_Vector_Utils.py

@@ -800,13 +800,11 @@ def ensure_MultiLineString(input_geometry):
     elif input_geometry.geom_type in ["LineString", "LinearRing"]:
         return geometry.MultiLineString([input_geometry])
     elif "Collection" in input_geometry.geom_type:
-        return geometry.MultiLineString(
-            (
-                line
-                for line in input_geometry.geoms
-                if line.geom_type in ["LineString", "LinearRing"]
-            )
-        )
+        valid_lines = []
+        for line in input_geometry.geoms:
+            if line.geom_type in ["LineString", "LinearRing"]:
+                valid_lines.append(line)
+        return geometry.MultiLineString(valid_lines)
     else:
         return geometry.MultiLineString()