0

I got an elevation dataset in .asc, I wanted to display it into Google Earth Pro so I made a Python script to convert it into GeoTIFF.

To my demise importing all the GeoTIFF at once only displays the last one; I have to import them one by one.

I notice that when Importing a GeoTIFF, Google Earth Pro does projection transformation. If you look at this image you will notice that the image is not in the orientation of the original bounding box

doing some automated reprojection (plane to sphere?)

If I add them one by one they perfectly line up

imported one by one

Then I added code to my script to make a KMZ out of all those GeoTIFF so I could import them in one click, however I had to convert the coordinate from Lambert93 to WGS84, and by doing so Google Earth Pro won't do the projection correction anymore giving me this:

my KMZ

How can I compute this reprojection computation in order to make a working KMZ?

The dataset is taken from here :

https://geoservices.ign.fr/bdalti (it's BDALTIV2_2-0_25M_ASC_LAMB93-IGN69_D014_2022-12-21)

and my Python code is like this:

import os
import numpy as np
from osgeo import gdal, osr
import matplotlib.pyplot as plt
import simplekml
import zipfile
import pyproj

source_folder = r"C:\checkout2\common-scripts\brgm\geoportail_carte_relief\source_files\BDALTIV2_2-0_25M_ASC_LAMB93-IGN69_D014_2022-12-21\BDALTIV2_2-0_25M_ASC_LAMB93-IGN69_D014_2022-12-21\BDALTIV2\1_DONNEES_LIVRAISON_2023-03-00115\BDALTIV2_MNT_25M_ASC_LAMB93_IGN69_D014"
output_folder = os.path.join(source_folder, "geotiffs")

if not os.path.exists(output_folder):
    os.makedirs(output_folder)

# Define the Lambert93 coordinate system
lambert93 = osr.SpatialReference()
lambert93.ImportFromEPSG(2154)
transformer = pyproj.Transformer.from_crs("EPSG:2154", "EPSG:4326", always_xy=True)

def apply_colormap(data_array, global_min, global_max):
    # Normalize the data using the global min and max
    norm_data = (data_array - global_min) / (global_max - global_min)

    # Apply the terrain colormap and then convert to 8-bit RGB
    colored = plt.cm.terrain(norm_data)
    return (colored[:, :, :3] * 255).astype(np.uint8)

# Step 1: Determine global min and max elevation values
global_min = float('inf')
global_max = float('-inf')

for filename in os.listdir(source_folder):
    if filename.endswith(".asc"):
        input_path = os.path.join(source_folder, filename)
        input_dataset = gdal.Open(input_path, gdal.GA_ReadOnly)
        data_array = input_dataset.ReadAsArray()

        current_min = np.nanmin(data_array)
        current_max = np.nanmax(data_array)

        global_min = min(global_min, current_min)
        global_max = max(global_max, current_max)

        input_dataset = None

for filename in os.listdir(source_folder):
    if filename.endswith(".asc"):
        input_path = os.path.join(source_folder, filename)
        output_path = os.path.join(output_folder, filename.replace(".asc", ".tif"))

        # Open the .asc file
        input_dataset = gdal.Open(input_path, gdal.GA_ReadOnly)

        # Get georeferencing information
        geotransform = input_dataset.GetGeoTransform()

        # Read the data and apply colormap
        data_array = input_dataset.ReadAsArray()
        rgb_array = apply_colormap(data_array, global_min, global_max)

        # Create an RGB output dataset
        driver = gdal.GetDriverByName('GTiff')
        output_dataset = driver.Create(output_path, input_dataset.RasterXSize, input_dataset.RasterYSize, 3,
                                       gdal.GDT_Byte)

        # Set georeferencing information
        output_dataset.SetGeoTransform(geotransform)
        output_dataset.SetProjection(lambert93.ExportToWkt())

        # Write the RGB data to the output dataset
        for band in range(3):
            output_dataset.GetRasterBand(band + 1).WriteArray(rgb_array[:, :, band])

        # Close datasets
        input_dataset = None
        output_dataset = None

# Create KMZ
kml = simplekml.Kml()

for filename in os.listdir(output_folder):
    if filename.endswith(".tif"):
        geotiff_path = os.path.join(output_folder, filename)

        dataset = gdal.Open(geotiff_path)
        geotransform = dataset.GetGeoTransform()

        xllcorner = geotransform[0]
        yllcorner = geotransform[3]
        xurcorner = xllcorner + geotransform[1] * dataset.RasterXSize
        yurcorner = yllcorner + geotransform[5] * dataset.RasterYSize

        xll_wgs84, yll_wgs84 = transformer.transform(xllcorner, yllcorner)
        xur_wgs84, yur_wgs84 = transformer.transform(xurcorner, yurcorner)

        groundoverlay = kml.newgroundoverlay(name=f"Elevation Overlay - {filename}")
        groundoverlay.icon.href = filename
        groundoverlay.latlonbox.north = yur_wgs84
        groundoverlay.latlonbox.south = yll_wgs84
        groundoverlay.latlonbox.east = xur_wgs84
        groundoverlay.latlonbox.west = xll_wgs84

temp_kml_path = os.path.join(output_folder, "combined_elevation.kml")
kml.save(temp_kml_path)

with zipfile.ZipFile(os.path.join(output_folder, "combined_elevation.kmz"), 'w') as kmz:
    kmz.write(temp_kml_path, "combined_elevation.kml")
    for geotiff_file in os.listdir(output_folder):
        if geotiff_file.endswith(".tif"):
            kmz.write(os.path.join(output_folder, geotiff_file), geotiff_file)

os.remove(temp_kml_path)

1 Answer 1

0

According to ChatGPT Lambert93 is for a flat surface and WGS84 for a sphere, using GDAL transformation did the trick:

for geotif_filename in geotif_filenames:
    # Reproject from Lambert-93 to WGS84
    lambert93 = osr.SpatialReference()
    lambert93.ImportFromEPSG(2154)

    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(4326)

    reprojection_tmp_file_path = "" # empty to use the MEM driver and not create a file on disc
    ds = gdal.Warp(reprojection_tmp_file_path, geotif_filename, format="MEM", srcSRS=lambert93, dstSRS=wgs84)

    # Convert reprojected GeoTIFF to PNG
    tmp_png = os.path.splitext(os.path.basename(geotif_filename))[0] + "_overlay.png"
    options = gdal.TranslateOptions(format='PNG', noData=0)
    gdal.Translate(tmp_png, ds, options=options)
    png_files.append(tmp_png)

    # Extract georeferencing info for KML
    geotransform = ds.GetGeoTransform()
    ulx, xres, xskew, uly, yskew, yres = geotransform
    lrx = ulx + (ds.RasterXSize * xres)
    lry = uly + (ds.RasterYSize * yres)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.