unable to create correct projection (TIF to PMTILES)

[idouiri9806]

Hi, I am currently using this code to generate PMTiles from TIF files. However, the projection appears zoomed in and significantly off. I'm not sure what's missing. Could you please take a look? This is a simplified version of what I am trying to achieve, as I want to apply color coding to terrain, etc . Also I already have a XYZ tile folder with png files, can i convert them to pmtiles directly ?

import os
import io
import rasterio
from rasterio.enums import Resampling
from rasterio.windows import Window
from PIL import Image
import subprocess
from multiprocessing import Pool, cpu_count
from pmtiles.writer import write
from pmtiles.tile import zxy_to_tileid, TileType, Compression

# Define default increments (in meters, approximately 100 feet)
default_green_increment = 30.48
default_yellow_increment = 30.48

def reproject_to_wgs84(input_path, output_path):
    """Reproject GeoTIFF to WGS84/EPSG:4326."""
    subprocess.run(['gdalwarp', '-t_srs', 'EPSG:4326', input_path, output_path])
    print(f"Reprojected {input_path} to {output_path}")

def convert_geotiff_to_pmtiles(geotiff_path, pmtiles_path, max_zoom=3):
    """Convert GeoTIFF to PMTiles."""
    print(f"Converting GeoTIFF {geotiff_path} to PMTiles {pmtiles_path}")
    try:
        with rasterio.open(geotiff_path) as src:
            with write(pmtiles_path) as writer:
                for z in range(max_zoom + 1):
                    tile_size = 256
                    num_tiles = 2 ** z
                    for x in range(num_tiles):
                        for y in range(num_tiles):
                            window = Window(x * tile_size, y * tile_size, tile_size, tile_size)
                            elevation_chunk = src.read(1, window=window, resampling=Resampling.bilinear)
                            img = Image.fromarray(elevation_chunk, 'L')
                            with io.BytesIO() as output:
                                img.save(output, format="PNG")
                                tileid = zxy_to_tileid(z, x, y)
                                writer.write_tile(tileid, output.getvalue())
                pmtiles_header = {
                    "min_zoom": 0,
                    "max_zoom": max_zoom,
                    "tile_type": TileType.PNG,
                    "tile_compression": Compression.NONE,
                }
                writer.finalize(pmtiles_header, {})
        print(f"PMTiles created at {pmtiles_path}")
    except Exception as e:
        print(f"Exception occurred while converting GeoTIFF to PMTiles: {e}")

def process_altitude(altitude, reprojected_geotiff_path, pmtiles_dir):
    pmtiles_output_path = os.path.join(pmtiles_dir, f"terrain_{altitude}.pmtiles")

    # Check if the output file already exists
    if os.path.exists(pmtiles_output_path):
        print(f"PMTiles for altitude {altitude} already exists. Skipping...")
        return

    # Convert the GeoTIFF to PMTiles
    convert_geotiff_to_pmtiles(reprojected_geotiff_path, pmtiles_output_path, max_zoom=3)

def worker(altitude):
    process_altitude(altitude, reprojected_geotiff_path, pmtiles_dir)

# Call the functions
geotiff_path = r"C:\Users\pannam\Desktop\pilotage_chart\world_terrain\output\us.tif"
reprojected_geotiff_path = r"C:\Users\pannam\Desktop\pilotage_chart\world_terrain\output\us_wgs84.tif"
output_dir = r"C:\Users\pannam\Desktop\pilotage_chart\world_terrain\output"

# Reproject the GeoTIFF to WGS84/EPSG:4326
reproject_to_wgs84(geotiff_path, reprojected_geotiff_path)

# Create directory for PMTiles
pmtiles_dir = os.path.join(output_dir, "us_terrain_pmtiles")
os.makedirs(pmtiles_dir, exist_ok=True)

# Determine the highest altitude from the GeoTIFF data
with rasterio.open(reprojected_geotiff_path) as src:
    max_altitude = src.read(1, masked=True).max()

# Generate PMTiles for different altitude thresholds with an increment of 30.48 meters (approximately 100 feet)
altitudes = range(0, 500, int(30.48))

# Use multiprocessing to process altitudes in parallel
if __name__ == "__main__":
    with Pool(processes=cpu_count()) as pool:
        pool.map(worker, altitudes)

[bdon]

Where are you converting to web mercator EPSG:3857?

Also I already have a XYZ tile folder with png files, can i convert them to pmtiles directly ?

https://github.com/protomaps/PMTiles/blob/main/python/bin/pmtiles-convert#L51