3

I create with C# a series of lat/lon value pair in the code(using i, j loop) and I want to extract the pixel value based on either lat/lon pair or pixel location (i, j).

I know there is one gdal utility called gdallocationinfo which can read pixel value very fast. But it is written in python.

I am wondering whether there is a way to use this gdal utility in my C# code, or other methods that can extract pixel value in a faster manner.

Final solution:

static void Main(string[] args)
    {
        Console.WriteLine(DateTime.Now.ToString());
        Gdal.AllRegister(); //Register all gdal data drives
        List<int> resolutionId = new List<int>() { 9 };

        Dataset ds = Gdal.Open("D:/SKnoZero.tif", Access.GA_ReadOnly);   //Read raster
        double[] gt = new double[6];
        ds.GetGeoTransform(gt); //Read geo transform info into array
        int Rows = ds.RasterYSize;  //Get the number of rows 
        int Cols = ds.RasterXSize;  //Get the number of columns
        Band band = ds.GetRasterBand(1);    //Read band 1
        double startX = gt[0];  //Upper left lon
        double startY = gt[3];  //Upper left lat
        double interval = gt[1];    //Cell size
        double x, y;    //Current lon and lat

        var filePath = "D:/SKcsv.csv";
        using (StreamWriter csv = new StreamWriter(filePath))
        {
            for (int k = 0; k < Rows; k++)  //read one line at a time
            {
                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++)  
                {
                    if (buf[r] != -32768)   //if pixel value is not NoData value
                    {
                        x = startX + r * interval;  //current lon                             
                        var geoidpair = GetVrgId(x, y, resolutionId);   //get geoid
                        var temp = geoidpair;
                        var newLine = string.Format("{0},{1},{2},{3}{4}", x, y, buf[r], temp[0], Environment.NewLine);
                        csv.Write(newLine);           
                     }
                 }
         }
         Console.WriteLine(DateTime.Now.ToString());
         Console.WriteLine("Press any key");
         Console.ReadLine();
    }
  • Have you successfully implemented GDAL in C#? – Michael Stimson Jul 31 '14 at 22:46
  • @MichaelMiles-Stimson Not yet. Other parts of my code have no relation with gdal. – shirleydd Jul 31 '14 at 22:51
  • 1
    gdallocationinfo is written in C not python. It uses RasterIO to read the pixel (line 475), this is available in C#. – Michael Stimson Jul 31 '14 at 22:51
  • @MichaelMiles-Stimson Oh gosh! This saves my life! Thank you! I only use it in OSGEO shell. Thought it is written in python. Is there source code in C# version? – shirleydd Jul 31 '14 at 22:58
  • The C# is a p-invoke layer written using SWIG. It works great. I use it in a significant application. – Reilly Aug 1 '14 at 10:43
3

I'm using OSGeo.GDAL for this code. It is extracted from my own working code but is incomplete, it will show you the basics of opening and reading a raster using GDAL.

The source for GDAL is downloadable, there is mostly C and C++ but some python... no C# that I know of - C# is windows only, GDAL is multi-platform. Many times I have referred to the source code for examples of what works and what doesn't. With C# you're mostly on your own as no really good documentation exists you can look at this though.

The ReadRaster parameters are StartRow, StartCol, ColumnsToRead, RowsToRead, BufferToStoreInto, BufferColumns, BufferRows, 0, 0. I hope that helps.

using OSGeo.GDAL; // way up the top

OSGeo.GDAL.Dataset pInDS;
pInDS = Gdal.Open(InImage.Text, Access.GA_ReadOnly);

// if the requested pixel is out of range return an error
// for that you need the range.
Rows = pInDS.RasterYSize;
Cols = pInDS.RasterXSize;

double[] pGT = new double[6];
pInDS.GetGeoTransform(pGT);
int col = (int) (X - pGT[0]) / pGT[1]
int row = (int) (pGT[3] - Y) / -pGT[5]

// pReadBlock is an array of the raster data type, you will need to duplicate that
// code dependent on each data type
OSGeo.GDAL.Band pInBand = pInDS.GetRasterBand(1);
bool IsThereAnError = pInBand.ReadRaster(row, col, 1, 1, pReadBlock, 1, 1, 0, 0);
Value = pReadBlock[0];

You might need to modify your reference path to include the GDAL bin/lib directory so it can find required libraries an I think you may need to link a specific lib GDAL_i.lib.

The geotransform array goes like this:

GeoTransform[0] = Xmin;             // Upper Left X
GeoTransform[1] = CellSize;         // W-E pixel size
GeoTransform[2] = 0;                // Rotation, 0 if 'North Up'
GeoTransform[3] = Ymax;             // Upper Left Y
GeoTransform[4] = 0;                // Rotation, 0 if 'North Up'
GeoTransform[5] = -CellSize;        // N-S pixel size

Note that the geotransform is from the upper left and cells are negative in the Y dimenstion. It is explained here.

  • This is awesome! Thank you Michael! I modified my code according to yours. I am wondering is there a way to check whether the raster is successfully read. When transferring the line starting with bool-----, there is an error saying that the type is not convertible between left and right sides. Because the value (-32768 as nodata value) I got from the code is not the true pixel value. I guess the raster is not successfully read? – shirleydd Jul 31 '14 at 23:35
  • My code: int[] buf = new int[1]; Band band = ds.GetRasterBand(1); band.ReadRaster(0, 0, ds.RasterXSize, ds.RasterYSize, buf, 1, 1, 0, 0); int value = buf[0]; Console.WriteLine(value); The value I got is wrong. – shirleydd Jul 31 '14 at 23:36
  • It's band.ReadRaster(row,col,1,1... you will need the geotransform to resolve X and Y to row and col, I will put that in. – Michael Stimson Jul 31 '14 at 23:41
  • 1
    @shirleydd: RTB. The C# API is nearly exactly the C++ API, so whatever you read on the gdal site applies. ReadRaster returns the same CPLErr as the C++. These are enums in C#. – Reilly Aug 1 '14 at 10:53
  • 1
    Why not read the raster one line (row) at a time and iterate through the values returned, instead of 1,1 for a pixel make it the raster cell width. GDAL also has a ReadBlock which is faster; some rasters have a block that is a row, others like tiled TIFF have blocks that are 128x128 but if you're reading the whole raster that shouldn't matter. – Michael Stimson Aug 4 '14 at 21:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.