0

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

  • It's not enough to state that coordinates are Cartesian -- you must specify which projected coordinate system (PCS) is in use (in ArcGIS, by defining the projection). Generally speaking, you should avoid on-the-fly projection of a raster by using its coordinate reference as the mapping canvas coordinate reference. – Vince Apr 11 '16 at 11:32
  • What country is your data supposed to be in? – Mintx Apr 11 '16 at 13:50
  • the country that this specific data set lies in is africa, inside the nothern region of kenya – jikoriko Apr 11 '16 at 15:52
  • Your code uses the Earth's radius, which is pretty inaccurate compared to what cs2cs can calculate. – AndreJ Apr 12 '16 at 5:56
1

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.

  • Thanks for the response, however the issue is not that I cant display the raster correctly inside QGIS, it does display correctly with the generated CRS, the issue I might not have been too clear on, is that I have my own custom map application, independent of QGIS, what I need is to be able to find out the wgs84 coordinate(latlon degrees) to be able to display a model I generate from the ascii grid data correctly inside the world. Using the method you mentioned could I extract the latitude/longitude after projection and replace the coordinates in the ascii file to the newly projected ones? – jikoriko Apr 11 '16 at 12:37
  • @jikoriko Please update your question to make this detail clear – Vince Apr 11 '16 at 14:18
  • With gdaltransform or cs2cs you can calculate the coordinates of the corner points from transverse mercator to degrees. But keep in mind that a rectangle in a projected CRS is not a rectangle in degrees anymore. And the cellsize in meters can not be recalculated to a fixed degree value. – AndreJ Apr 11 '16 at 16:09
  • You can get the projected coordinates of the other corners from xllcorner, yllcorner, cellsize, nrows and ncols, then perform the same transformation on them. Note that xllcorner/yllcorner might be mid-cell, or edge of the cell. – AndreJ Apr 12 '16 at 6:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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