[OSGeo-UK] Fw: Get Lat/Long/Altitude data from GeoTiff File
Paul Harwood
runette at gmail.com
Sat Feb 20 13:46:09 PST 2021
I may be able to help - I have done a fair amount of this in c#.
First a general note - the c# API is very close to the Java API - so in
general you can look at that documentation.
Next - kind of a criticism. Your code is kind of messy and uncommented -
for instance, you seem to get the GeoTransform twice for each row - which
is massively not optimized (just once for the raster will do) - and you get
the value from gt before you get it and I don't think that will work. If
you want help from a community list then it is a good idea to have clear
and well-commented code.
As best as I can tell - you want to convert each value in a geotiff band
into x,y,z coordinates in BNG coordinates. If that is not a correct
understanding then much of the rest of this mail will be wrong.
You can only get z from a geotiff if it is a DEM or DTM - but I assume that
you know what you are doing there. You do have a problem with projection
that I will come back to.
Finally - what you are doing is not easy - especially with a reprojection.
Where I do this same thing in C# I actually use a different library (PDAL -
see https://pdal.io/ using the pdalc library with extensions into c#
https://github.com/PDAL/CAPI). This uses GDAL to load the GeoTIFF but can
reproject and return the data as a cloud of points in an array. If you want
to know about that - let me know and I will pm you.
Basically - you want it in this order :
1 Get the Geotransform (once)
2. Project the GeoTransform to BNG
3 iterate through the raster
1 - Get the Geotransform - you know how to do that
2 - Reproject
A couple of points:
- It is better to import the destination CRS from an EPSG code than from a
proj string - there are some gaps in proj strings and they are being phased
out. To get BNG as the destination - I would do ;
SpatialReference dst = new SpatialReference(null);
dst.importFromEPSG(27700);
- you have not mentioned the CRS of the geotiff. I assume that it is
correctly defined in the GeoTIFF o GDAL should pick it up. I also
assume that it is EPSG:4326 (i.e. WGS84 latlong). In your sample code, you
seem to have the src and dst the wrong way around for what you said in the
text.
Create the top left and bottom right corners from the GT. these will be
4326 - convert them as you have done but you only need to this for the two
corners - not each point. Note that the raster is always 2d.
The pixel sizes from the gt are in degrees (because 4326). You probably
want to dump them and recreate them. 27700 (i.e. bng) is projected and use
1m units - so you can just do a simple interpolation from the corners.
3 Iterate
At this point, you can get the 2D coordinates for each point by LERPing
from the corners. You asked about the z coord.
You have to aware that neither GDAL nor GeoTIFF 'know' about the 3rd
dimension. To them it i just data. If the tiff is a DEM or DTm, then the
band value for the pixel is the z coord.
However, since it is just data it is not reprojected in any way. so what
you get will depend on the pc of the geotiff.
For instance, if the tiff is shipped EPSG:4326 or wgs84 then the values
will often be an orthogonal height above the WG84 ellipsoid. After this
process, they still will be. However, BNG does not use the WGS84 ellipsoid
- it uses OSGB36. So - a naive user might expect the z values to be
orthogonal height above OSGB36 or above OSGM15 (i.e. UK MSW datum). You
just need to b aware.
On Fri, 19 Feb 2021 at 15:50, Wael El-Sayegh <wael_elsayegh at hotmail.com>
wrote:
>
>
>
>
>
> I am trying to obtain and convert latitude and longitude and altitude
> values into British National Grid coordinates from GeoTIFF using GDAL C#
> file, I managed to read the coordinate system and corner coordinates as
>
> The thing is, since that GDAL is more of Python documentation and I can't
> go through it, I am not sure what is the next step or how to do it using
> SpatialReference() and CoordinateTransformation().
>
> This is the code I tried so far but still getting the wrong numbers, I
> think its because they are projected coordinate system and needed to be
> converted in British coordinate system.
>
> How can I do this?
>
>
>
> PS: I am not sure if the XY coordinates are right, also not sure how to
> read Z coordinates.
>
>
>
> Looking forward to hearing from you
>
>
>
> GdalConfiguration.ConfigureGdal();
>
> GdalConfiguration.ConfigureOgr();
>
> Gdal.AllRegister();
>
>
>
> Dataset ds = Gdal.Open(fileName, Access.GA_ReadOnly);
>
>
>
>
>
>
>
>
>
> double[] gt = new double[6];
>
> int Rows = ds.RasterYSize;
>
> int Cols = ds.RasterXSize;
>
> double startX = gt[0];
>
> double startY = gt[3];
>
> double interval;
>
> Band band = ds.GetRasterBand(1);
>
> double xx, y;
>
>
>
> for (int k = 0; k < Rows; k++)
>
> {
>
>
>
> ds.GetGeoTransform(gt);
>
>
>
> interval = gt[1];
>
> ds.GetGeoTransform(gt);
>
> y = startY - k * interval; //Current lat
>
> int[] buf = new int[Cols];
>
> //ReadRaster parameters are StartCol, StartRow,
> ColumnsToRead, RowsToRead, BufferToStoreInto, BufferColumns, BufferRows, 0,
> 0
>
> band.ReadRaster(0, k, Cols, 1, buf, Cols, 1, 0, 0);
>
> //iterate each item in one line
>
> for (int r = 0; r < Cols; r++)
>
> {
>
> row = XYZData.NewRow();
>
> if (buf[r] != -32768) //if pixel value is not NoData
> value
>
> {
>
> xx = startX + r * interval; //current lon
>
>
>
> SpatialReference srs = new
> SpatialReference("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601271625
> +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs");
>
>
>
> SpatialReference dst = new SpatialReference("");
>
> dst.ImportFromProj4("+proj=longlat +ellps=WGS84
> +datum=WGS84 +no_defs");
>
>
>
>
>
> double[] p = new double[3];
>
> p[0] = xx;
>
> p[1] = y;
>
> p[2] = 0;
>
>
>
> var coordinateTransform = new
> CoordinateTransformation(srs, dst);
>
>
>
> coordinateTransform.TransformPoint(p);
>
>
>
> var latLong = new LatitudeLongitude(p[1], p[0]);
>
> var cartesian = GeoUK.Convert.ToCartesian(new
> Wgs84(), latLong);
>
>
>
> row["X"] = cartesian.X;
>
> row["Y"] = cartesian.Y;
>
> row["Z"] = 0;
>
> XYZData.Rows.Add(row.ItemArray);
>
>
>
>
>
> }
>
>
>
> _______________________________________________
> UK mailing list
> UK at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/uk
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/uk/attachments/20210220/a65e4d41/attachment-0001.html>
More information about the UK
mailing list