2

I am trying to use libtiff.net to read elevation data from a GeoTIFF file. So far I have mostly just been able to read metadata from the file using the example at libtiff.net's webpage.

But howto read elevation data I do not understand... I tried first reading with Tiff.ReadScanline() as described here but the file I have seems to be stored differently (probably in tiles if I understand it correctly)

Here is the metadata (as far as I have been able to read) (the tiff file is from the danish terrain elevation data set):

Tiff c:\Users***\DTM_1km_6170_500.tif, page 0 has following tags set:

IMAGEWIDTH System.Int32 : 2500

IMAGELENGTH System.Int32 : 2500

BITSPERSAMPLE System.Int16 : 32

COMPRESSION BitMiracle.LibTiff.Classic.Compression : ADOBE_DEFLATE

PHOTOMETRIC BitMiracle.LibTiff.Classic.Photometric : MINISBLACK

STRIPOFFSETS System.UInt64[] : System.UInt64[]

SAMPLESPERPIXEL System.Int16 : 1

STRIPBYTECOUNTS System.UInt64[] : System.UInt64[]

PLANARCONFIG BitMiracle.LibTiff.Classic.PlanarConfig : CONTIG

PREDICTOR BitMiracle.LibTiff.Classic.Predictor : FLOATINGPOINT

TILEWIDTH System.Int32 : 256

TILELENGTH System.Int32 : 256

TILEOFFSETS System.UInt64[] : System.UInt64[]

TILEBYTECOUNTS System.UInt64[] : System.UInt64[]

SAMPLEFORMAT BitMiracle.LibTiff.Classic.SampleFormat : IEEEFP

DATATYPE System.Int16 : 3

GEOTIFF_MODELPIXELSCALETAG System.Int32 : 3 GEOTIFF_MODELPIXELSCALETAG System.Byte[] : Ù?Ù?

GEOTIFF_MODELTIEPOINTTAG System.Int32 : 6 GEOTIFF_MODELTIEPOINTTAG System.Byte[] : A ^WA

34735 System.Int32 : 36 34735 System.Byte[] : ± ± #° èd )#

34736 System.Int32 : 3 34736 System.Byte[] :

34737 System.Int32 : 30 34737 System.Byte[] : ETRS89 / UTM zone 32N|ETRS89|

42113 System.Int32 : 6 42113 System.Byte[] : -9999

The code I have written so far is as follows:

namespace GeoTIFFReader
{
  public class GeoTIFF
  {
    private double[,] heightmap;
    private double dx;
    private double dy;
    private double startx;
    private double starty;


    public GeoTIFF(string fn)
    {
      using (Tiff tiff = Tiff.Open(fn, "r"))
      {
        if (tiff == null)
        {
          // Error - could not open
          return;
        }



        int width = tiff.GetField(TiffTag.IMAGEWIDTH)[0].ToInt();
        int height = tiff.GetField(TiffTag.IMAGELENGTH)[0].ToInt();
        heightmap = new double[width, height];
        FieldValue[] modelPixelScaleTag = tiff.GetField(TiffTag.GEOTIFF_MODELPIXELSCALETAG);
        FieldValue[] modelTiePointTag = tiff.GetField(TiffTag.GEOTIFF_MODELTIEPOINTTAG);

        byte[] modelPixelScale = modelPixelScaleTag[1].GetBytes();
        dx = BitConverter.ToDouble(modelPixelScale, 0);
        dy = BitConverter.ToDouble(modelPixelScale, 8) * -1;

        byte[] modelTransformation = modelTiePointTag[1].GetBytes();
        double originLon = BitConverter.ToDouble(modelTransformation, 24);
        double originLat = BitConverter.ToDouble(modelTransformation, 32);

        startx = originLon + dx / 2.0;
        starty = originLat + dy / 2.0;

        double curx = startx;
        double cury = starty;

        FieldValue[] bitsPerSampleTag = tiff.GetField(TiffTag.BITSPERSAMPLE);

        FieldValue[] tilewtag = tiff.GetField(TiffTag.TILEWIDTH);
        FieldValue[] tilehtag = tiff.GetField(TiffTag.TILELENGTH);
        int tilew = tilewtag[0].ToInt();
        int tileh = tilehtag[0].ToInt();

        var tile = new byte[tilew*tileh];

        //var scanline = new byte[tiff.ScanlineSize()]; Does not work... wrong format
        for (int il = 0; il < height; il++)
        {
          //tiff.ReadScanline(scanline, il); // Load il'th line of data 
          for (int ir = 0; ir < width; ir++)
          {

            // Here I would like to read each pixel data that contains elevation in gray-scale in f32 as far I as I understand from metadata

            //object value = scanline[ir];
            //heightmap[ir, il] = double.Parse(value.ToString());
          }
        }

        Console.WriteLine(heightmap.ToString());
      }

    }
  }
}

So if anyone knows howto extract this data, that would be very much appreciated.

Erik Thysell
  • 1,160
  • 1
  • 14
  • 31

4 Answers4

2

So I stumbled across some hinting that lead me to find an answer to the specific question..:

    int tileSize = tiff.TileSize();
    for (int iw = 0; iw < nWidth; iw += tilew)
    {
      for (int ih = 0; ih < nHeight; ih += tileh)
      {
        byte[] buffer = new byte[tileSize];
        tiff.ReadTile(buffer, 0, iw, ih, 0, 0);
        for (int itw = 0; itw < tilew; itw++)
        {
          int iwhm = ih + itw;
          if (iwhm > nWidth - 1)
          {
            break;
          }
          for (int ith = 0; ith < tileh; ith++)
          {
            int iyhm = iw + ith;
            if (iyhm > nHeight - 1)
            {
              break;
            }
            heightMap[iwhm, iyhm] =
              BitConverter.ToSingle(buffer, (itw * tileh + ith) * 4);
          }
        }
      }
    }

EDIT 2018-09-20: @Graviton - Sorry for the long response time.. but here is the complete class I have used with a constructor that takes the filename as input and populates the heightMap (but @Nazonokaizijin looks a bit nicer and slimmer):

using System;
using BitMiracle.LibTiff.Classic;
using System.IO;

namespace GeoTIFFReader
{
  public class GeoTIFF
  {
    private float[,] heightMap;

    public float[,] HeightMap
    {
      get { return heightMap; }
      private set { heightMap = value; }
    }
    private int nWidth;

    public int NWidth
    {
      get { return nWidth; }
      private set { nWidth = value; }
    }
    private int nHeight;

    public int NHeight
    {
      get { return nHeight; }
      private set { nHeight = value; }
    }
    private double dW;

    public double DW
    {
      get { return dW; }
      private set { dW = value; }
    }
    private double dH;

    public double DH
    {
      get { return dH; }
      private set { dH = value; }
    }
    private double startW;

    public double StartW
    {
      get { return startW; }
      private set { startW = value; }
    }
    private double startH;

    public double StartH
    {
      get { return startH; }
      private set { startH = value; }
    }


    public GeoTIFF(string fn)
    {
      using (Tiff tiff = Tiff.Open(fn, "r"))
      {
        if (tiff == null)
        {
          // Error - could not open
          return;
        }

        nWidth = tiff.GetField(TiffTag.IMAGEWIDTH)[0].ToInt();
        nHeight = tiff.GetField(TiffTag.IMAGELENGTH)[0].ToInt();
        heightMap = new float[nWidth, nHeight];
        FieldValue[] modelPixelScaleTag = tiff.GetField(TiffTag.GEOTIFF_MODELPIXELSCALETAG);
        FieldValue[] modelTiePointTag = tiff.GetField(TiffTag.GEOTIFF_MODELTIEPOINTTAG);

        byte[] modelPixelScale = modelPixelScaleTag[1].GetBytes();
        dW = BitConverter.ToDouble(modelPixelScale, 0);
        dH = BitConverter.ToDouble(modelPixelScale, 8) * -1;

        byte[] modelTransformation = modelTiePointTag[1].GetBytes();
        double originLon = BitConverter.ToDouble(modelTransformation, 24);
        double originLat = BitConverter.ToDouble(modelTransformation, 32);

        startW = originLon + dW / 2.0;
        startH = originLat + dH / 2.0;

        FieldValue[] tileByteCountsTag = tiff.GetField(TiffTag.TILEBYTECOUNTS);
        long[] tileByteCounts = tileByteCountsTag[0].TolongArray();

        FieldValue[] bitsPerSampleTag = tiff.GetField(TiffTag.BITSPERSAMPLE);
        int bytesPerSample = bitsPerSampleTag[0].ToInt() / 8;


        FieldValue[] tilewtag = tiff.GetField(TiffTag.TILEWIDTH);
        FieldValue[] tilehtag = tiff.GetField(TiffTag.TILELENGTH);
        int tilew = tilewtag[0].ToInt();
        int tileh = tilehtag[0].ToInt();

        int tileWidthCount = nWidth / tilew;
        int remainingWidth = nWidth - tileWidthCount * tilew;
        if (remainingWidth > 0)
        {
          tileWidthCount++;
        }

        int tileHeightCount = nHeight / tileh;
        int remainingHeight = nHeight - tileHeightCount * tileh;
        if (remainingHeight > 0)
        {
          tileHeightCount++;
        }

        int tileSize = tiff.TileSize();
        for (int iw = 0; iw < nWidth; iw += tilew)
        {
          for (int ih = 0; ih < nHeight; ih += tileh)
          {
            byte[] buffer = new byte[tileSize];
            tiff.ReadTile(buffer, 0, iw, ih, 0, 0);
            for (int itw = 0; itw < tilew; itw++)
            {
              int iwhm = ih + itw;
              if (iwhm > nWidth - 1)
              {
                break;
              }
              for (int ith = 0; ith < tileh; ith++)
              {
                int iyhm = iw + ith;
                if (iyhm > nHeight - 1)
                {
                  break;
                }
                heightMap[iwhm, iyhm] =
                  BitConverter.ToSingle(buffer, (itw * tileh + ith) * 4);
              }
            }
          }
        }
      }
    }
  }
}
Erik Thysell
  • 1,160
  • 1
  • 14
  • 31
  • I tried to get your code running, but unable because the code is not strictly compilable. Would you want to make your answer compilable and self-contained? That would help a lot – Graviton Jun 19 '18 at 08:54
1

You can do it something like this:

private void ReadTiff()
{
    Tiff tiff = Tiff.Open("myfile.tif", "r");
    if (tiff == null)
        return;

    //Get the image size
    int imageWidth = tiff.GetField(TiffTag.IMAGEWIDTH)[0].ToInt();
    int imageHeight = tiff.GetField(TiffTag.IMAGELENGTH)[0].ToInt();

    //Get the tile size
    int tileWidth = tiff.GetField(TiffTag.TILEWIDTH)[0].ToInt();
    int tileHeight = tiff.GetField(TiffTag.TILELENGTH)[0].ToInt();
    int tileSize = tiff.TileSize();

    //Pixel depth
    int depth = tileSize / (tileWidth * tileHeight);

    byte[] buffer = new byte[tileSize];

    for (int y = 0; y < imageHeight; y += tileHeight)
    {
        for (int x = 0; x < imageWidth; x += tileWidth)
        {
            //Read the value and store to the buffer
            tiff.ReadTile(buffer, 0, x, y, 0, 0);

            for (int kx = 0; kx < tileWidth; kx++)
            {
                for (int ky = 0; ky < tileHeight; ky++)
                {
                    //Calculate the index in the buffer
                    int startIndex = (kx + tileWidth * ky) * depth;
                    if (startIndex >= buffer.Length)
                        continue;

                    //Calculate pixel index
                    int pixelX = x + kx;
                    int pixelY = y + ky;
                    if (pixelX >= imageWidth || pixelY >= imageHeight)
                        continue;

                    //Get the value for the target pixel
                    double value = BitConverter.ToSingle(buffer, startIndex);
                }
            }
        }
    }

    tiff.Close();
}
Nazonokaijin
  • 121
  • 10
  • Welcome to Stack Overflow! Please don't answer just with source code. Try to provide a nice description about how your solution works. See: [How do I write a good answer?](https://stackoverflow.com/help/how-to-answer). Thanks – sɐunıɔןɐqɐp Sep 20 '18 at 06:20
  • I've downloaded multiple geoTIFF files - to make sure the files themselves weren't the problem - and every file I run using your code, or similar code crashes whenever the values from TiffTag.TILEWIDTH or TiffTag.TILELENGTH are used. Any idea what's causing these values to be null when reading the geoTIFF files? – Jimmy Nov 01 '18 at 18:54
  • Some tiff files don't have the Tile structure. In your case, it seems that the tiff files you downloaded don't have the tile structure. (This is the reason why your code crashed.) You will need to get the pixel data using different approach than above. – Nazonokaijin Nov 15 '18 at 05:25
  • I am getting error in line double value = BitConverter.ToSingle(innerBuffer, startIndex); Destination array is not long enough... (startIndex = 131070, buffer.Length = 131072). So I added if (startIndex > innerBuffer.Length - 4) – Martin86 Jul 13 '21 at 08:46
1

It has a lib (in c#) complementary to LibTiff that makes the elevation query given a latitude/longitude, available at GeoTiffCOG:

  GeoTiff geoTiff = new GeoTiff(file_tiff);
  double value = geoTiff.GetElevationAtLatLon(latitude, longitude);

It also supports Cloud Optimized GeoTiff (COG) by adding URI as an argument.

Fabric.io
  • 11
  • 1
0

I believe the problem is the PREDICTOR. Instead of placing the data in the file, it is LZW encoding of the differences in the data, with FLOATINGPOINT predictor. This was proprietary with Adobe. I have been searching for code do decrypt PREDICTOR_FLOATINGPOINT, myself. I think this github link may have library code that will read and decode it.

Fabulous
  • 2,393
  • 2
  • 20
  • 27