[GIS] Finding WGS84 coordinate from ASCII grid and calculating scale over a slippy tile map

coordinate systemqgis

I am currently trying to calculate the WGS84 coordinate from an ESRI ASCII grid file, which the xllCorner and yllCorner coordinates are I believe in Cartesian coordinates, here is a sample of the data's header:

ncols         3744
nrows         2110
xllcorner     5000160.447
yllcorner     -5304137.100
cellsize      90
nodata_value  -9999

the aim of converting the coordinate is so that I can display a model generated from the file at the correct world position and correct scale inside my own custom map application, which uses a slippy tile system.

I have used QGIS to sample the data as a raster and it displayed correctly using this generated CRS:

+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601272 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs

I don't know if it is possible to use the data visible in that CRS to get the values that I need, but from using the obvious values like lat_0 and lon_0 I can't calculate the correct location, I also noticed that the projection says Transverse Mercator, but plugging in the raw values to a UTM calculator did not work.

Edit:

I have played around with getting a WGS84 coordinate from AndreJ's suggestion by re-projecting the data into the correct CRS, then looking at the layers properties inside QGIS I have managed to find a reasonable coordinate for the data's origin, however the issue that I now have is with the scaling of the model to fit correctly over my slippy tile system, here are some screenshots to show the issue:

without water plane

with water plane

just to clarify a few things:

  • the height data is being composed into a 3D mesh, where the z value is the height value
  • I render the mesh fully transparent, then render a water plane at a given z/height value, the visible water pixels are when the z values exceed the z value of the meshes transparent pixels (using the depth buffer in OpenGL)
  • the mesh stores its vertices x meters apart from each other, where the number of meters is the cellsize given by the ASCII grid.

I set the water plane to be at a reasonable height around the same height the visible lake is at, as you can see the outline of the lake does not match up with the mesh, due to scale.

Here is some of the code I am using to calculate positions and scales:

public static GpsCoordinate OffsetByMeters(GpsCoordinate start, double metersNorth, double metersEast)
    {
        GpsCoordinate result = new GpsCoordinate(start.Latitude, start.Longitude);

        double metersLat = metersNorth / EarthRadius;
        double metersLon = metersEast / (EarthRadius * Math.Cos(Math.PI * start.Latitude / 180));
        result.Latitude += metersLat * 180 / Math.PI;
        result.Longitude += metersLon * 180 / Math.PI;

        return result;
    }

public static Vector3 LatLonToFloatXy(double latitude, double longitude, int zoom)
    {
        latitude = Wrap(latitude, MinLatitude, MaxLatitude);
        longitude = Wrap(longitude, MinLongitude, MaxLongitude);

        double x = (longitude + 180) / 360;
        double sinLatitude = Math.Sin(latitude * Math.PI / 180);
        double y = 0.5 - Math.Log((1 + sinLatitude) / (1 - sinLatitude)) / (4 * Math.PI);

        uint mapSize = MapSize(zoom);
        float floatX = (float)Wrap(x * mapSize + 0.5, 0, mapSize - 1);
        float floatY = (float)Wrap(y * mapSize + 0.5, 0, mapSize - 1);

        Vector3 position = new Vector3(floatX, floatY, 0);
        return position;
    }

void Render()
{
    Vector3 start = LatLonToFloatXy(_startCoord.Latitude, _startCoord.Longitude, CURRENTMAPZOOM);
    GpsCoordinate endCoord = OffsetByMeters(_startCoord, -_meshHeight, _meshWidth);
    //mesh width and height are in meters
    Vector3 end = LatLonToFloatXy(endCoord.Latitude, endCoord.Longitude, CURRENTMAPZOOM);
    float scaleFactor = (end - start).X / meshWidth;
    Vector3 scale = Vector3.One * scaleFactor;
    //now I Render the mesh and plane with the start position and scale
}

So to summarise I am calculating the scale by getting an offset GPS coordinate in meters by the size of mesh in meters, converting the start GPS coordinate and end GPS coordinate to pixel coordinates, and measuring the difference in them. Using this method I hope to get a correctly scaled height model over the map imagery, but there must be an issue with my implementation somewhere.

Also I have ended up editing the ASCII grid file to contain the WGS84 latitude and longitude, however it would be more desirable to be able to calculate the coordinates to make the import of data easier in future. I am not sure as to how easy this would be as I assume the ASCII grid files coordinates can be in many CRS's

Best Answer

Since you have the CRS parameters, you can assign that CRS to your data with gdal_translate:

gdal_translate -a_srs "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601272 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +unity=m +no_defs" src_dataset temp.tif

and then use gdalwarp to reproject the data to latlon degrees:

gdalwarp -t_srs EPSG:4326 temp.tif out.tif

UTM is a transverse mercator projection, with a false Easting of 500km to a numbered 6-degree longitude raster. Your custom CRS does not fit that rule.

Once you have assigned the CRS, QGIS should be able to show the raster on the right spot against any basemap.

Related Question