Date   
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).

Re: Rasterio result different than gdal_calc

anand@...
 

Hi Sean,

Sure, I will try those changes you have mentioned. here is the gdal_calc expression I used :

gdal_calc.py -A T1_B5.TIF --A_band=1 -B T1_B4.TIF --B_band=1 --outfile=ndvi_gal_1.TIF --calc='(A.astype(float)-B)/(A.astype(float)+B)' --type='Float32'
--co='COMPRESS=LZW'

Re: Rasterio result different than gdal_calc

Sean Gillies
 

Hi,

Can you also include the expressions that you're using with gdal_calc? Without them we can only guess at the cause of the difference in the results. I don't see any big problems with the way you're using Rasterio. I have a couple comments below:

On Wed, Nov 21, 2018 at 3:35 PM anand via Groups.Io <anand=geospoc.com@groups.io> wrote:
I am trying to recreate gdal_calc for simple calculation using below code.

import numpy
import rasterio
def merge_bands(calculationType, band1Image, band2Image, outputFilePath):
# Read raster bands directly to Numpy arrays.
#
# We handle the connections with "with"
with rasterio.open(band1Image) as dsRed:
bandRed = dsRed.read(1)

You could convert this to a float32 now instead of doing it twice in a statement below.

 
with rasterio.open(band2Image) as dsNIR:
bandNIR = dsNIR.read(1)

Same comment.

 
# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')

Unless you're doing this in gdal_calc too, this could be a cause of differences. I think it would be better to add keyword argument `masked=True` when you read the source bands. This will give you masked arrays that are safer to work when calculating the NDVI.

ndvi = (bandNIR.astype(float)-bandRed.astype(float))/(bandNIR.astype(float) + bandRed.astype(float))

If you didn't specify `astype(float)` here, Numpy would convert to float64 as soon as it performed the division operation, and would return a float64 ndarray.
 
kwargs = dsRed.meta
kwargs.update(dtype=rasterio.float32, count=1, compress='lzw')
with rasterio.open(outputFilePath, 'w', **kwargs) as dst:
dst.write_band(1, ndvi.astype(rasterio.float32))

But the result of this script is different than gdal_calc:

gdal_cal result range -0.1 - +0.09
rasterio result range 0.02 - 0.1

Like I said, we need to see the gdal_calc expression. I suspect nodata pixels are involved.

--
Sean Gillies

Rasterio result different than gdal_calc

anand@...
 

I am trying to recreate gdal_calc for simple calculation using below code.

import numpy
import rasterio
def merge_bands(calculationType, band1Image, band2Image, outputFilePath):
# Read raster bands directly to Numpy arrays.
#
# We handle the connections with "with"
with rasterio.open(band1Image) as dsRed:
bandRed = dsRed.read(1)
 
with rasterio.open(band2Image) as dsNIR:
bandNIR = dsNIR.read(1)
 
# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')
ndvi = (bandNIR.astype(float)-bandRed.astype(float))/(bandNIR.astype(float) + bandRed.astype(float))
kwargs = dsRed.meta
kwargs.update(dtype=rasterio.float32, count=1, compress='lzw')
with rasterio.open(outputFilePath, 'w', **kwargs) as dst:
dst.write_band(1, ndvi.astype(rasterio.float32))

But the result of this script is different than gdal_calc:

gdal_cal result range -0.1 - +0.09
rasterio result range 0.02 - 0.1
 

Re: Unclear how to implement image alignment based on GCPs

Sean Gillies
 

Hi,

On Tue, Nov 6, 2018 at 10:04 PM <mmweene15@...> wrote:
Hi, thanks for your response. Let me provide some background. I have two sets of images that capture the same scene but were taken using separate cameras. This causes some degree of GPS discrepancy between the two sets of images, and the task is to align them based on GCPs generated either manually or automatically (irrelevant for this problem). Maybe I have completely misunderstood the purpose of the `reproject` function. I have included a sample of my GCPs below:

What an interesting problem. The reproject() function is mainly for cartographic reprojection, but it might be useful in your case.
 

`GroundControlPoint(row=array(357), col=array(682), x=array(22.67637871), y=array(-33.94166626), z=0.0, id='45e0b74c-54a4-48c8-b876-1dfcc184c6c8')`

The row and column reference the position in the aligned image that corresponds to the x and y coordinates in the misaligned image. Perhaps my understanding of how to construct this is wrong as well? The documentation doesn't seem to cover this case all that clearly.

The row, col, x, and y for a GroundControlPoint need to be floats, not arrays of floats. I'm sorry that the documentation of this class is inadequate!


Thanks for your help.

You're welcome! I think that adjusting the GCPs may help, but may not be entirely enough.

--
Sean Gillies

Re: Bug: 'JPEG2000' is not a valid Compression

Sean Gillies
 

Hi, 

This is indeed a bug. I would be much obliged if you reported it in the Rasterio issue tracker.

I haven't used NITF in a long time. I see in https://www.gdal.org/frmt_nitf.html that there's a lack of symmetry in the compression creation option name ("IC"), values ("C3", "M3", "C8"), and what gets stored in the metadata ("COMPRESSION" and "JPEG2000"). This will complicate things a bit, more than in the GeoTIFF situation where I mostly work.


On Wed, Nov 7, 2018 at 12:17 PM <nicholas.maxwell@...> wrote:
Is this the right place to report a bug? I can open an issue on github if you like.

I'm reading a NITF, and getting this error. Should add 'JPEG2000' to the Compression enum?

In[14]: src.tags(ns='IMAGE_STRUCTURE')
Out[14]: {'COMPRESSION': 'JPEG2000'}

In[15]: src.compression
Traceback (most recent call last):
  File "rasterio/_base.pyx", line 840, in rasterio._base.DatasetBase.compression.__get__
  File "/lego/lib/python3.7/enum.py", line 307, in __call__
    return cls.__new__(cls, value)
  File "/lego/lib/python3.7/enum.py", line 555, in __new__
    return cls._missing_(value)
  File "/lego/lib/python3.7/enum.py", line 568, in _missing_
    raise ValueError("%r is not a valid %s" % (value, cls.__name__))
ValueError: 'JPEG2000' is not a valid Compression

thanks



--
Sean Gillies

Bug: 'JPEG2000' is not a valid Compression

nicholas.maxwell@...
 

Is this the right place to report a bug? I can open an issue on github if you like.

I'm reading a NITF, and getting this error. Should add 'JPEG2000' to the Compression enum?

In[14]: src.tags(ns='IMAGE_STRUCTURE')
Out[14]: {'COMPRESSION': 'JPEG2000'}

In[15]: src.compression
Traceback (most recent call last):
  File "rasterio/_base.pyx", line 840, in rasterio._base.DatasetBase.compression.__get__
  File "/lego/lib/python3.7/enum.py", line 307, in __call__
    return cls.__new__(cls, value)
  File "/lego/lib/python3.7/enum.py", line 555, in __new__
    return cls._missing_(value)
  File "/lego/lib/python3.7/enum.py", line 568, in _missing_
    raise ValueError("%r is not a valid %s" % (value, cls.__name__))
ValueError: 'JPEG2000' is not a valid Compression

thanks

Re: Make _boundless_vrt_doc public

Guy Doulberg
 

Hi Sean

I can't find the dev-group to discuss this issue, where can I find it?

Thanks, Guy

Re: Unclear how to implement image alignment based on GCPs

mmweene15@...
 

Hi, thanks for your response. Let me provide some background. I have two sets of images that capture the same scene but were taken using separate cameras. This causes some degree of GPS discrepancy between the two sets of images, and the task is to align them based on GCPs generated either manually or automatically (irrelevant for this problem). Maybe I have completely misunderstood the purpose of the `reproject` function. I have included a sample of my GCPs below:

`GroundControlPoint(row=array(357), col=array(682), x=array(22.67637871), y=array(-33.94166626), z=0.0, id='45e0b74c-54a4-48c8-b876-1dfcc184c6c8')`

The row and column reference the position in the aligned image that corresponds to the x and y coordinates in the misaligned image. Perhaps my understanding of how to construct this is wrong as well? The documentation doesn't seem to cover this case all that clearly.

Thanks for your help.

Re: Unclear how to implement image alignment based on GCPs

Sean Gillies
 

Hi,

I'm having trouble understanding the intent of code. It looks like you are calling reproject() with source and destination arrays of the same shape, with the same coordinate reference system on each side of the reprojection. I'm not what we should expect in this situation; we are not currently testing for this. Can you write back and show your GCPs?

Unclear how to implement image alignment based on GCPs

mmweene15@...
 

I am trying to replicated the gdal functionality for performing image alignment based on GCPs. The equivalent gdal commands are the following:

`gdal_translate -of GTiff -gcp GCP1 -gcp GCP2... src_raster dst_raster`
`gdalwarp -s_srs S_SRS -t_srs T_SRS -tps dst_raster dst_raster2`

This would create a new raster, `dst_raster2` that has been reprojected according to the specified GCPs. It is unclear, based on the documentation, how to achieve this using the rasterio.warp module. I have tried using the rasterio.warp.reproject function to no avail. I've included a snippet of code below. The resulting image does not match the image produced by gdal. Any guidance would be appreciated. Note that the variables are descriptively named.

with rasterio.open(target_image_path, 'r') as f:
target_profile = f.profile
target_image = image_tools.rasterio_read_bands(f)
print(target_profile)

destination = np.empty(target_image.shape, dtype=target_profile['dtype'])
resolution = (target_profile['transform'].a, target_profile['transform'].e)
transform, width, height = rasterio_warp.calculate_default_transform(src_crs,
dst_crs,
gcps=gcps,
resolution=resolution,
height=target_profile['height'],
width=target_profile['width'])

rasterio_warp.reproject(source=target_image,
destination=destination,
gcps=gcps,
src_crs=src_crs,
dst_crs=src_crs,
# src_transform=target_profile['transform'],
dst_transform=target_profile['transform'],
resampling=getattr(enums.Resampling, resampling_method),
num_threads=num_threads)

Re: Make _boundless_vrt_doc public

Sean Gillies
 

Hi Denis,

I'm not prepared to make it public yet. I'm not confident that the signature of the method is finished and I believe we should have a discussion on the dev group about VRTs and how much support we'll have for them in Rasterio's public API.

For now, I recommend copying it into your own project.


On Wed, Oct 31, 2018 at 9:39 AM Denis Rykov <rykovd@...> wrote:
Hi,

We are going to use rasterio.vrt._boundless_vrt_doc for making a XML VRT out of dataset in our library.
But we are little afraid because the name of this function tells that it is a private one.
Are there any reasons for this? What do you think about to make this function public?



--
Sean Gillies

Make _boundless_vrt_doc public

Denis Rykov
 

Hi,

We are going to use rasterio.vrt._boundless_vrt_doc for making a XML VRT out of dataset in our library.
But we are little afraid because the name of this function tells that it is a private one.
Are there any reasons for this? What do you think about to make this function public?

Re: Any way to get/set the raster band scale and offset?

Kris Vanhoof
 

That would be great, thanks a lot!

Re: Any way to get/set the raster band scale and offset?

Sean Gillies
 

Hi Kris,

Raster scale and offset don't surface in Rasterio versions 1.0.9 and older only because no user has spoken up them. I'd be happy to add these for 1.0.10.


On Mon, Oct 29, 2018 at 4:22 PM Kris Vanhoof <kris.vanhoof@...> wrote:
Hi! I am using rasterio to convert a number of geotiffs that were originally
created using GDAL.  To preserve the image metadata I'm currently using a combination of
the 'profile', 'descriptions' and 'units' properties of DatasetBase as well as 'tags()'.

However, it is unclear to me how to get and set the 'Scale' and 'Offset' value for a raster band.

Since there is a 'units' property in 'DatasetBase', I would have expected there to be a 'scales' and 'offsets'.

Likewise, in 'DatasetWriterBase' there is a 'set_band_unit()', but no 'set_band_offset()' or 'set_band_scale()'.

The GDAL functions that I would normally use are:

GDALGetRasterUnitType()
GDALGetRasterScale()
GDALGetRasterOffset()

GDALSetRasterUnitType()
GDALSetRasterScale()
GDALSetRasterOffset()

Is there any way to access these values using rasterio?



--
Sean Gillies

Any way to get/set the raster band scale and offset?

Kris Vanhoof
 

Hi! I am using rasterio to convert a number of geotiffs that were originally
created using GDAL.  To preserve the image metadata I'm currently using a combination of
the 'profile', 'descriptions' and 'units' properties of DatasetBase as well as 'tags()'.

However, it is unclear to me how to get and set the 'Scale' and 'Offset' value for a raster band.

Since there is a 'units' property in 'DatasetBase', I would have expected there to be a 'scales' and 'offsets'.

Likewise, in 'DatasetWriterBase' there is a 'set_band_unit()', but no 'set_band_offset()' or 'set_band_scale()'.

The GDAL functions that I would normally use are:

GDALGetRasterUnitType()
GDALGetRasterScale()
GDALGetRasterOffset()

GDALSetRasterUnitType()
GDALSetRasterScale()
GDALSetRasterOffset()

Is there any way to access these values using rasterio?

Rasterio 1.0.9

Sean Gillies
 

Hi all.

Rasterio 1.0.9 is on PyPI now. It has a bunch of bug fixes that are especially important if you are deploying systems on AWS or using the WarpedVRT class.


Thanks for the bug reports, analysis, and patience, everybody!

--
Sean Gillies

Re: Windowed reads and coordinates

Koshy Thomas
 

Turns out the transform object needs to actually be

transform = rasterio.windows.transform(window, dataset.transform)

for windowed reads. 

Windowed reads and coordinates

Koshy Thomas
 

Hi all

I'm trying to create vector features out of raster files using the `features.shapes` module and I'm running into a weird issue where when I run a windowed read on a `.tif`, the Y coordinates of the polygons are thrown off while the X coordinates seem correct. If I run the same logic without a windowed read (against the whole tif), it seems to work correctly. Any ideas? I'm attaching the .tif and the code screenshot is below:

Re: GDAL version policy

Denis Rykov
 

Thank you Sean!

I made a PR at rasterio-wheels repository which updates GDAL to the latest stable version.

Besides this we discovered that we need a rasterio wheels built against a master version of GDAL. What do you advise us in this case? Does it make sense having such wheels on PyPI?