Fw: Get Lat/Long/Altitude data from GeoTiff File

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Fw: Get Lat/Long/Altitude data from GeoTiff File

Wael El-Sayegh



 

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
[hidden email]
https://lists.osgeo.org/mailman/listinfo/uk
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Get Lat/Long/Altitude data from GeoTiff File

Yakus
You're probably going to get more assistance with this if you post it on the gdal-dev list
http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html

Thanks

Matt

On Fri, 19 Feb 2021 at 15:50, Wael El-Sayegh <[hidden email]> 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
[hidden email]
https://lists.osgeo.org/mailman/listinfo/uk

_______________________________________________
UK mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/uk
Reply | Threaded
Open this post in threaded view
|

Re: Fw: Get Lat/Long/Altitude data from GeoTiff File

Paul Harwood
In reply to this post by Wael El-Sayegh
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 :

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 <[hidden email]> 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
[hidden email]
https://lists.osgeo.org/mailman/listinfo/uk

_______________________________________________
UK mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/uk