[gdal-dev] C# bindings for (fast) reproject, resample, and send to DB

[mladen-g at distributel.net]

Hello,

I've been using GDAL for a while, but I am still a newbie to the API.   
I see there are some experts on this group, so apologies if I'm asking  
something simple (my searches did not reveal the answers I was looking  
for).

Here is my problem:

I need to take an input raster, split into areas of various sizes,  
resample each area into a 256x256 tile, reproject to EPSG:4326 and  
send the resulting pixels to a DB.

GDAL is perfect for this, but I am using the C# bindings and I'm not  
sure what is the most performant way to do it.

I used the following two references to create something that works,  
but I am wondering if it can be improved:
https://gis.stackexchange.com/questions/297831/using-gdalwarp-with-c-bindings
https://lists.osgeo.org/pipermail/gdal-dev/2017-February/046046.html

How I do it right now:

1)  Copy the input grid to RAM disk for fast access
2)  Invoke Gdal.wrapper_GDALWarpDestName() for reprojection,  
resampling, and output to a /vsimem/ memory-mapped file
3)  Read the memory-mapped file into a Dataset (returned by  
wrapper_GDALWarpDestName())
4)  Write the Dataset values to DB

What I would like to know:

1)  To skip the /vsimem file creation, how can I use  
Gdal.wrapper_GDALWarpDestDs() instead of -DestName?  It requires a  
dataset as a parameter, but I am not sure how to initialize it  
correctly.  I've tried Gdal.GetDriverByName("VRT").Create(...) and  
.GetDriverByName("MEM").Create(...), but all the values end up being 0  
when I pass the dataset to wrapper_GDALWarpDestDs()?

2)  Is there a more direct, efficient way to do what I'm doing?  Here  
is the relevant code

<snip>
using (Dataset rasterInput = Gdal.Open(inputRasterPath, Access.GA_ReadOnly))
{
   IntPtr[] ptr = {Dataset.getCPtr(inputRasterPath).Handle};
   GCHandle gcHandle = GCHandle.Alloc(ptr, GCHandleType.Pinned);
...
    var dss = new  
SWIGTYPE_p_p_GDALDatasetShadow(gcHandle.AddrOfPinnedObject(), false,  
null);
..
    foreach (var tileExtent in overlappingTileExtents) // tileExtent  
just contains minX/minY/maxX/maxY
    {
    ...

       string vsiMemFilePath =  
$"/vsimem/{destinationTableName}#{extentName}#tile{indexNumber}.vrt";
       string[] options = {
                                 "-of", "VRT",
                                 "-srcnodata", "99999.9",
                                 "-dstnodata", "99999.9",
                                 "-ot", "Float32",
                                 "-t_srs", "EPSG:4326",
                                 "-r", "average"),
                                 "-te", tileExtent.xmin.ToString(),  
tileExtent.ymin.ToString(), tileExtent.xmax.ToString(),  
tileExtent.ymax.ToString(),
                                 "-ts", "256", "256"
                             };

        using (Dataset ds =  
Gdal.wrapper_GDALWarpDestName(vsiMemFilePath, 1, dss, new  
GDALWarpAppOptions(options), null, null))
        {
...
          //read the geotransform + first band and write each pixel  
lat/long + value to the DB
...
        }

     Gdal.Unlink(vsiMemFilePath)

....

Regards,
Mladen