What are the corresponding fields to set in libtiff or libgeotiff, given a minimal ESRI ASCII file?

Question

I have an ESRI ASCII file of the form:

ncols         5 
nrows         4 
xllcorner     0 
yllcorner     0 
cellsize      10 
NODATA_value  -9999 
25.4 26.1 27 28.6 27.7
25 26 26.4 27.9 27.4
25.1 25.8 26.8 28.6 27.6
27.5 28 27.7 30.6 28.3

And I need to use libtiff.net ( or its C++ equivalent, libtiff, or libgeotiff, or GDAL, or any other C# or C++ libraries) to convert it to a geotiff file.

But I simply don't know what are the fields to set, there are a lot of fields such as TIFFTAG_IMAGEWIDTH, TIFFTAG_SAMPLESPERPIXEL or TIFFTAG_BITSPERSAMPLE which I don't know whether they are relevant or not.

Given an ESRI ASCII file as per above, how to use the libtiff.net or libtiff library to create a geotiff file?

Why am I doing this?

In my application, I have a mesh, and I would have to convert the mesh to a raster file. I want to create such a geotiff raster file directly from the mesh sampling, without creating an intermediate ASCII file first because the intermediate file size is very big compared to the final geotiff output file.

The only approach I can think of is to directly manipulate the geotiff file using libtiff.net while I am obtaining the grid points on the mesh.

How to do this, or is there a better approach?

Answer 1

Even if you are planning not to use GDAL in your application (more on that later), nothing prevents you from generating a geotiff from the ESRI ASCII GRID (gdal_translate -of "GTiff" in.asc out.tif), and inspect the TIFF tags of the generated file (these are the necessary tags to generate a geotiff from the given grid).

AsTiffTagViewer gives the following output (similar to TiffTags utility output):

TagCode (Count DataType): Value  // my comments

ImageWidth (1 Short): 5    // ncols
ImageLength (1 Short): 4   // nrows
BitsPerSample (1 Short): 32
Compression (1 Short): Uncompressed
Photometric (1 Short): MinIsBlack
StripOffsets (1 Long): 260
SamplesPerPixel (1 Short): 1
RowsPerStrip (1 Short): 4    //nrows
StripByteCounts (1 Long): 80
PlanarConfig (1 Short): Contig
SampleFormat (1 Short): 3
33550 (3 Double): 
33922 (6 Double): 
42113 (6 ASCII): -9999    // NODATA_value

As you can see there are 11 standard Tiff tags and 3 non standard tags (but we know the Data Types and more importantly the dimensions of that last 3 tags, 3,6,6). Let's confirm we have 2 GeoTiff tags and one non standard GDAL-specific tag.

Libgeotiff C library is distributed with listgeo utility for dumping GeoTIFF metadata. The output:

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                0                0                
         0                40               0                
      ModelPixelScaleTag (1,3):
         10               10               0                
      End_Of_Tags.
   Keyed_Information:
      End_Of_Keys.
   End_Of_Geotiff.

Corner Coordinates:
Upper Left    (       0.000,      40.000)
Lower Left    (       0.000,       0.000)
Upper Right   (      50.000,      40.000)
Lower Right   (      50.000,       0.000)
Center        (      25.000,      20.000)

By the dimensions of the Tagged_Information we can identify the following 2 tags. Also, since the grid is regular (equal X and Y spacing, and there are no skewed grid lines) we can establish the following formulas:

33550 tag:

ModelPixelScaleTag = [ cellsize , cellsize , 0 ]

33922 tag:

ModelTiepointTag = [ 0 , 0 , 0 , UpperLeftCorner_X , UpperLeftCorner_Y , 0] 

    where
            UpperLeftCorner_X = xllcorner 

            UpperLeftCorner_Y = yllcorner + cellsize * nrows

That leaves the last tag 42113. The geotiff format does not have a standard tag for a nodata value. GDAL stores band nodata value in the non standard TIFFTAG_GDAL_NODATA ASCII tag (code 42113).

Finally, as an example, we can write a function that relates the header of the grid (ncols, nrows, cellsize, xllcorner, yllcorner) with the Tiff tags, using Libgeotiff C library:

void SetUpTIFFDirectory(TIFF *tif)
{
    double tiepoints[6];
    double pixscale[3];

    double upperLeftCorner_X, upperLeftCorner_Y;

    upperLeftCorner_X = xllcorner;

    upperLeftCorner_Y = yllcorner + (cellsize*nrows);

    tiepoint[0] = 0.0;
    tiepoint[1] = 0.0;
    tiepoint[2] = 0.0;
    tiepoint[3] = upperLeftCorner_X;
    tiepoint[4] = upperLeftCorner_Y;
    tiepoint[5] = 0.0;

    pixscale[0] = cellsize;
    pixscale[1] = cellsize;
    pixscale[2] = 0.0;

    TIFFSetField(tif,TIFFTAG_IMAGEWIDTH,    ncols);
    TIFFSetField(tif,TIFFTAG_IMAGELENGTH,   nrows);
    TIFFSetField(tif,TIFFTAG_BITSPERSAMPLE, 32);
    TIFFSetField(tif,TIFFTAG_COMPRESSION,   COMPRESSION_NONE);
    TIFFSetField(tif,TIFFTAG_PHOTOMETRIC,   PHOTOMETRIC_MINISBLACK);
    TIFFSetField(tif,TIFFTAG_STRIPOFFSETS,  260L);
    TIFFSetField(tif,TIFFTAG_SAMPLESPERPIXEL, 1);
    TIFFSetField(tif,TIFFTAG_ROWSPERSTRIP, nrows;
    TIFFSetField(tif,TIFFTAG_STRIPBYTECOUNTS, 80L);
    TIFFSetField(tif,TIFFTAG_PLANARCONFIG,  PLANARCONFIG_CONTIG);
    TIFFSetField(tif,TIFFTAG_SAMPLEFORMAT,  3;

    TIFFSetField(tif,GTIFF_TIEPOINTS, 6,tiepoints);
    TIFFSetField(tif,GTIFF_PIXELSCALE, 3,pixscale);

}

Note: When you say you can't use GDAL, there's an In Memory Raster format that can probably be used as a temporary placeholder for you raster while you are adding mesh samples:

https://www.gdal.org/frmt_mem.html

Answer 2

Actually this task can be easily done by using GDAL library ( or its .Net equivalent, GDAL.Net). There is even an example here in Python:

ncols         174
nrows         115
xllcorner     14.97
yllcorner     -34.54
cellsize      0.11

And the Python script:

if __name__ == '__main__':
# Import libs
import numpy, os
from osgeo import osr, gdal

# Set file vars
output_file = "out.tif"

# Create gtif
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, 174, 115, 1, gdal.GDT_Byte )
raster = numpy.zeros( (174, 115) )

# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform( [ 14.97, 0.11, 0, -34.54, 0, 0.11 ] )

# set the reference info 
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84")
dst_ds.SetProjection( srs.ExportToWkt() )

# write the band
dst_ds.GetRasterBand(1).WriteArray(raster)

Converting the code to .Net should be trivial enough.