GeoTIFF libtiff.net 在 C# 中获取高程数据

GeoTIFF libtiff.net get elevation data in c#

我正在尝试使用 libtiff.net 从 GeoTIFF 文件中读取高程数据。 到目前为止,我基本上只能使用 libtiff.net 网页上的示例从文件中读取元数据。

但是我不明白如何读取高程数据...我尝试按照 here 所述首先使用 Tiff.ReadScanline() 进行读取,但我的文件似乎存储方式不同(如果我理解正确)

这是元数据(据我所知)(tiff 文件来自丹麦地形高程数据集):

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

目前我写的代码如下:

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());
      }

    }
  }
}

因此,如果有人知道如何提取这些数据,将不胜感激。

所以我偶然发现了一些提示,这些提示让我找到了特定问题的答案..:[=​​13=]

    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);
          }
        }
      }
    }

编辑 2018-09-20: @Graviton - 很抱歉响应时间长..但这是完整的 class 我使用了一个构造函数,该构造函数将文件名作为输入并填充 heightMap (但@Nazonokaizijin 看起来更好一点并且苗条):

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);
              }
            }
          }
        }
      }
    }
  }
}

我认为问题出在 PREDICTOR。它不是将数据放在文件中,而是使用 FLOATINGPOINT 预测器对数据差异进行 LZW 编码。这是 Adob​​e 专有的。我自己一直在寻找解密 PREDICTOR_FLOATINGPOINT 的代码。我觉得this github link 可能有可以读取和解码它的库代码。

你可以这样做:

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();
}

它有一个库(在 c# 中)与 LibTiff 互补,可以在给定 latitude/longitude 的情况下进行海拔查询,可在 GeoTiffCOG:

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

它还通过添加 URI 作为参数来支持云优化 GeoTiff (COG)。