Topics

GDALRasterizeGeometries - Function Signature

jdcasta@...
 

Hi this is a little bit  of a strange question about how rasterio is using the GDALRasterizeGeometries function. This c-function is called in _features.pyx here:

with InMemoryRaster(image=image, transform=transform) as mem:
  exc_wrap_int(
  GDALRasterizeGeometries(
  mem.handle(), 1, mem.band_ids,num_geoms, geoms, NULL,
  mem.gdal_transform, pixel_values, options, NULL, NULL))

My question has to do with the bolded arguments, 6 & 7. The python function docstring expresses the notion that the affine matrix, transform, is used convert between the shapes coordinate system to the pixels provided in image. However looking at the CPP source code for this function I see this function declaration:


CPLErr GDALRasterizeGeometries ( GDALDatasetH  hDS,
    int  nBandCount,
    int *  panBandList,
    int  nGeomCount,
    OGRGeometryH *  pahGeometries,
    GDALTransformerFunc  pfnTransformer,
    void *  pTransformArg,
    double *  padfGeomBurnValue,
    char **  papszOptions,
    GDALProgressFunc  pfnProgress,
    void *  pProgressArg 
  )

As you can see its it seems rasterio is passing NULL to pfnTransformer and the affine transform to pTransformArg. This pTransformArg is not really meant to be used for the affine matrix, but as an argument to a valid pfnTransformer.

It *think* everything still works because of this code here that checks the arguments. If it sees that pfnTransformer is NULL it will automatically try to get the transform from the raster dataset. But I have not idead how in memory dataset will have that information
so I keep thinking this theory must be wrong.
if( pfnTransformer == nullptr )
  {
  bNeedToFreeTransformer = true;
   
  char** papszTransformerOptions = nullptr;
  double adfGeoTransform[6] = { 0.0 };
  if( poDS->GetGeoTransform( adfGeoTransform ) != CE_None &&
  poDS->GetGCPCount() == 0 &&
  poDS->GetMetadata("RPC") == nullptr )
  {
  papszTransformerOptions = CSLSetNameValue(
  papszTransformerOptions, "DST_METHOD", "NO_GEOTRANSFORM");
  }
   
  pTransformArg =
  GDALCreateGenImgProjTransformer2( nullptr, hDS,
  papszTransformerOptions );
  CSLDestroy( papszTransformerOptions );
   
  pfnTransformer = GDALGenImgProjTransform;
  if( pTransformArg == nullptr )
  {
  return CE_Failure;
  }
  }


In the end I am just trying to figure out how in the world you are passing just the affine transform into the  GDALRasterizeGeometries as  pTransformArg.

PS. I'm used to github formatting so this has been a little different to type up!
In the end I was just trying to write some CPP code that is similar to rasterio for rasterizing freatures (in memory).

Sean Gillies
 

You're right, the transform that Rasterio passes to GDALRasterizeGeometries is not used. Worse, there looks like risk of corruption. I'm going to change Rasterio to pass NULL instead.

When GDALRasterizeGeometries is called without a transformer, the following code is executed: https://github.com/OSGeo/gdal/blob/master/gdal/alg/gdalrasterize.cpp#L686-L710. Therein you can see how a transformer is constructed and configured.


On Mon, Dec 3, 2018 at 8:52 PM <jdcasta@...> wrote:
Hi this is a little bit  of a strange question about how rasterio is using the GDALRasterizeGeometries function. This c-function is called in _features.pyx here:

with InMemoryRaster(image=image, transform=transform) as mem:
  exc_wrap_int(
  GDALRasterizeGeometries(
  mem.handle(), 1, mem.band_ids,num_geoms, geoms, NULL,
  mem.gdal_transform, pixel_values, options, NULL, NULL))

My question has to do with the bolded arguments, 6 & 7. The python function docstring expresses the notion that the affine matrix, transform, is used convert between the shapes coordinate system to the pixels provided in image. However looking at the CPP source code for this function I see this function declaration:


CPLErr GDALRasterizeGeometries ( GDALDatasetH  hDS,
    int  nBandCount,
    int *  panBandList,
    int  nGeomCount,
    OGRGeometryH *  pahGeometries,
    GDALTransformerFunc  pfnTransformer,
    void *  pTransformArg,
    double *  padfGeomBurnValue,
    char **  papszOptions,
    GDALProgressFunc  pfnProgress,
    void *  pProgressArg 
  )

As you can see its it seems rasterio is passing NULL to pfnTransformer and the affine transform to pTransformArg. This pTransformArg is not really meant to be used for the affine matrix, but as an argument to a valid pfnTransformer.

It *think* everything still works because of this code here that checks the arguments. If it sees that pfnTransformer is NULL it will automatically try to get the transform from the raster dataset. But I have not idead how in memory dataset will have that information
so I keep thinking this theory must be wrong.
if( pfnTransformer == nullptr )
  {
  bNeedToFreeTransformer = true;
   
  char** papszTransformerOptions = nullptr;
  double adfGeoTransform[6] = { 0.0 };
  if( poDS->GetGeoTransform( adfGeoTransform ) != CE_None &&
  poDS->GetGCPCount() == 0 &&
  poDS->GetMetadata("RPC") == nullptr )
  {
  papszTransformerOptions = CSLSetNameValue(
  papszTransformerOptions, "DST_METHOD", "NO_GEOTRANSFORM");
  }
   
  pTransformArg =
  GDALCreateGenImgProjTransformer2( nullptr, hDS,
  papszTransformerOptions );
  CSLDestroy( papszTransformerOptions );
   
  pfnTransformer = GDALGenImgProjTransform;
  if( pTransformArg == nullptr )
  {
  return CE_Failure;
  }
  }


In the end I am just trying to figure out how in the world you are passing just the affine transform into the  GDALRasterizeGeometries as  pTransformArg.

PS. I'm used to github formatting so this has been a little different to type up!
In the end I was just trying to write some CPP code that is similar to rasterio for rasterizing freatures (in memory).



--
Sean Gillies

jdcasta@...
 

Yeah I thought so. FYI the code you linked is the same one I embedded above. It looks like it really is just relying upon the in memory raster already having a configured transform associated with it.  I did the same to get it to work like so:

        GDALDatasetH memDS = GDALCreate( GDALGetDriverByName("MEM"), proc_num_str, cols, rows, 0, GDT_Byte, NULL );
        // Set transform
        reinterpret_cast<GDALDataset*>(memDS)->SetGeoTransform(transform);

Here are set the in memory raster dataset to use the specified transform (double [6]). Now I can call GDALRasterizeGeometries with NULL for pTransform and pTransformArg and it is handled appropriately.

Thanks!