Modifying many_points.rs example to render elevation points

[BogdanOlar]

I'm trying to render some elevation points, but I'm pretty sure I'm using the wrong coordinate system(s).

The code below replaces the original generate_points() in many_points.rs example, and it essentially stacks a bunch of points near the N00E00 coordinates.

fn generate_points() -> Vec<ColoredPoint> {
    let fpath = Path::new("galileo/examples/data/N45E024.hgt");
    // Hardcoded lat/lon to match the data in the file
    let lat = 45.0;
    let lon = 24.0;
    // Distance in degrees between each elevation point
    let res_inc = 1.0 / 3600.0;
    // File contains 3601*3601 elevation values, 2 bytes each (i16)
    const GRID_COUNT: usize = 3601;
    const ELEVATION_VALUE_BYTE_COUNT: usize = 2;

    let file = File::open(fpath).expect("Cannot read file.");
    let mut data = Vec::new();
    let byte_count = BufReader::new(file).read_to_end(&mut data).unwrap();
    // GRID_COUNT columns, GRID_COUNT lines, and each elevation value is 2 bytes (big endian)
    assert_eq!(
        byte_count,
        GRID_COUNT * GRID_COUNT * ELEVATION_VALUE_BYTE_COUNT
    );

    let points: Vec<ColoredPoint> = data
        .iter()
        .step_by(ELEVATION_VALUE_BYTE_COUNT)
        .enumerate()
        .zip(data[1..].iter().step_by(ELEVATION_VALUE_BYTE_COUNT))
        .map(|((i, b0), b1)| {
            let ei: i16 = (((*b0 as u16) << 8usize) | (*b1 as u16)) as i16;
            let xi = i % GRID_COUNT;
            let yi = i / GRID_COUNT;

            let x = lat + (res_inc * (xi as f64));
            let y = lon - (res_inc * (yi as f64));
            let z = ei as f64;

            // FIXME: this is definitely wrong
            // let x = x * 20_000.0;
            // let y = y * 20_000.0;

            let red = z * (255 as f64) / 2600.0;
            let red = red as u8;
            let color = Color::rgba(red, 0, 0, 200);

            ColoredPoint {
                point: Point3d::new(x, y, z),
                color,
            }
        })
        .step_by(10)
        .collect();

    let highest_point = points
        .iter()
        .max_by(|p0, p1| p0.point.z().total_cmp(&p1.point.z()))
        .unwrap();

    log::info!("Generated {} points", points.len());
    log::info!("First point: {}", &points[0].point);
    log::info!("Highest point: {}", highest_point.point);

    // [2024-09-01T16:54:31Z INFO  elevation_points] Generated 1296721 points
    // [2024-09-01T16:54:31Z INFO  elevation_points] First point: {45, 24, 330}
    // [2024-09-01T16:54:31Z INFO  elevation_points] Highest point: {45.736111111111114, 24.4, 2516}

    points
}

When uncommenting the FIXME lines above, I can (kind of) see the height map of the area I'm interested in, with the small inconvenience that it's placed off the coast of Africa, in the Atlantic ocean, instead of Romania.

Do I need to create GeoPoint2d points instead of Point3d, and then somehow convert them to Point2d and add another dimension to get a Point3d?

Not even sure this question makes sense, so any advice is much appreciated.

N45E024.hgt.zip

[Maximkaaa]

Hey, @BogdanOlar . At least you can see the mountains, that is already cool! :-)

To project your points into correct coordinate system:

1. Get the projection for Web Mercator CRS (the one that the map uses)

let projection = Crs::EPSG3857.get_projection();

2. For every point in your data set create a 2d geo point and project it to your map's CRS:

let lat_p = lat + (res_inc * (xi as f64));
let lon_p = lon - (res_inc * (yi as f64));
let point = GeoPoint2d::latlon(lat, lon);
let projected = projection.project(&point);

3. Create a 3d-point by adding elevation to the projected point.

let point3d = Point3d::new(projected.x, projected.y, z);

This will give you points in the right position of the world, but a couple things I would also check in your code:

1. You add x coordinate to the latitude. I'm not familiar with the hgt format, but x is usually associated with the longitude, not latitude. So make sure that rows in that data format represent latitude increments.
2. I would expect the elevation data to be in meters. But the coordinates in Mercator projection are actually pseudo-meters. Putting z coordinates as is to the map would make the mountains look flatter than they really are. You can easily adjust your height by multiplying height in meters by 1 / cos(lat).

Let me know if this helps.

[BogdanOlar]

Yes! The tile is rendered correctly now!

Here's the updated code, in case anyone else has the same problem:

fn generate_points() -> Vec<ColoredPoint> {
    // File contains 3601*3601 elevation values, 2 bytes each (i16)
    const GRID_COUNT: usize = 3601;
    const ELEVATION_VALUE_BYTE_COUNT: usize = 2;
    // GRID_COUNT columns, GRID_COUNT lines, and each elevation value is 2 bytes (big endian)
    const DATA: &[u8; GRID_COUNT * GRID_COUNT * ELEVATION_VALUE_BYTE_COUNT] =
        include_bytes!("data/N45E024.hgt");
    // Initial lat/lon are specified in the file name, and corresponds to the lower-left corner of the tile
    const INIT_LAT: f64 = 45.0;
    const INIT_LON: f64 = 24.0;
    const LAT_LON_RESOLUTION: f64 = 1.0 / 3600.0;
    const MAX_ELEVATION: f64 = 2600.0;

    let projection = Crs::EPSG3857
        .get_projection()
        .expect("Could not get projection");

    let points: Vec<ColoredPoint> = DATA
        .iter()
        .step_by(ELEVATION_VALUE_BYTE_COUNT)
        .zip(DATA[1..].iter().step_by(ELEVATION_VALUE_BYTE_COUNT))
        .enumerate()
        .step_by(10)
        .map(|(i, (b0, b1))| {
            let ei: i16 = (((*b0 as u16) << 8usize) | (*b1 as u16)) as i16;
            let xi = i % GRID_COUNT;
            let yi = GRID_COUNT - (i / GRID_COUNT);

            let lon = INIT_LON + ((xi as f64) * LAT_LON_RESOLUTION);
            let lat = INIT_LAT + ((yi as f64) * LAT_LON_RESOLUTION);
            let elevation = ei as f64;

            let geo_point: GeoPoint2d = NewGeoPoint::latlon(lat, lon);
            let projected: Point2d = projection.project(&geo_point).unwrap();
            let point = Point3d::new(projected.x, projected.y, elevation / lon.cos());

            let red = (elevation * (u8::MAX as f64) / MAX_ELEVATION) as u8;
            let color = Color::rgba(red, 10, 0, 200);

            ColoredPoint { point, color }
        })
        .collect();

    let highest_point = points
        .iter()
        .max_by(|p0, p1| p0.point.z().total_cmp(&p1.point.z()))
        .unwrap();

    log::info!("Generated {} points", points.len());
    log::info!("First point: {}", &points[0].point);
    log::info!("Highest point: {}", highest_point.point);

    // [2024-09-02T17:19:30Z INFO  elevation_points] Generated 1296721 points
    // [2024-09-02T17:19:30Z INFO  elevation_points] First point: {2671667.7790385657, 5780393.73445169, 777.9734364313445}
    // [2024-09-02T17:19:30Z INFO  elevation_points] Highest point: {2678068.649759179, 5703715.594948082, 4489.821225163708}

    points
}

Thanks so much! I'm learning a lot of new stuff here.

[BogdanOlar]

Also, from what I understand, hgt is not really a format.

The data was taken from the NASA Shuttle Radar Topography Mission (SRTM) dataset (it's available for free, but you have to register for an account).

Each file contains 3601*3601, 16-bit elevation values, expressed in meters, and the geographical coordinates of the tile are encoded in the filename.

So, yeah, not really a file format.

I had to figure out that the first 16-bit value you read does not correspond to the coordinates in the filename, but to the lat+1, lon coordinates. :)

Edit: more info here

[Maximkaaa]

Looks super cool!

[BogdanOlar]

I wonder how difficult it would be to take those PNG tiles received from open street map and mesh them over those points (is mesh the correct term?). I would think it's not trivial.

[pka]

That would be awesome and I think, it shouldn't be that difficult. Example from Maplibre: https://maplibre.org/maplibre-gl-js/docs/examples/3d-terrain/
A world-wide tileset is hosted on AWS: https://registry.opendata.aws/terrain-tiles/ with the URL-Template https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png
See also this georaster example for reading a DTM over HTTP, but georaster is currently only supporting TIFF formats.

[Maximkaaa]

Yes, it is not too difficult if we talk only about raster tiles. Vector tiles and feature layers bring much more serious challenge...