Switch to warped vrt internally for boundless reads #1161

sgillies · 2017-10-02T21:30:55Z

Rasterio has offered "boundless reads" as a feature – you can read an image window (in pixel space) that goes beyond the source raster's rows and columns and rasterio performs VRT-like source reading and compositing. The math has been buggy: at best it doesn't exactly match GDAL's math and, at worst, can drop a line or row of data.

This PR fixes the bug by switching from an approximation of GDAL's VRT logic to the use of a temporary VRT that represents the boundless window. It provides more accurate and reliable results and eliminates 2 chunks of brittle code in _io.pyx.

I'm using a warped VRT because GDAL doesn't surface a non-warped VRT constructor in its C API. As far as I can tell, reprojection is short-circuited when the input and output crs are the same. I think this is good enough for 1.0 and we could look at a more optimal implementation after that.

I'm concerned that source overviews will be skipped when we use the VRT like this. Notice that we read from the VRT at its full resolution. It seems that the test at https://github.com/OSGeo/gdal/blob/8372431040f2b17bb5f6b643c2bfdfcec8d6f630/gdal/frmts/vrt/vrtsourcedrasterband.cpp#L169 will thus fail and source overviews would be skipped. I'm going to ask GDAL developers for confirmation and guidance on this point.

sgillies · 2017-10-03T18:51:58Z

Note this is somewhat related to #1131, in which we've found the limits of very large warped VRTs.

rouault · 2017-10-03T19:58:48Z

Short answer: no, source overviews will not be used that way

Long answer: I see you use GDALCreateWarpedVRT(), and this function will not try to use sourcce overviews.

Preliminary steps: create dataset of size 2048x2048, with uniform grey overview of size 1024x1024

gdal_translate http://svn.osgeo.org/gdal/trunk/autotest/gcore/data/byte.tif test.tif -outsize 2048 2048 --config GDAL_SKIP DODS
gdal_translate test.tif test.tif.ovr -outsize 1024 1024 -scale 0 255 127 127

Now test.c that basically use GDALCreateWarpedVRT() to create a view of reduced resolution at 1024x1024 pixels

#include "gdal_vrt.h"
#include "gdalwarper.h"
#include "cpl_port.h"

int main()
{
    GDALAllRegister();
    GDALDatasetH hSrcDS = GDALOpen("test.tif", GA_ReadOnly);
    double adfGeoTransform[6];
    GDALGetGeoTransform(hSrcDS, adfGeoTransform);
    adfGeoTransform[1] *= 2;
    adfGeoTransform[5] *= 2;
    GDALWarpOptions *psWOptions = GDALCreateWarpOptions();
    psWOptions->hSrcDS = hSrcDS;
    psWOptions->nBandCount = 1;
    psWOptions->panSrcBands = (int*)CPLMalloc(1*sizeof(int));
    psWOptions->panDstBands = (int*)CPLMalloc(1*sizeof(int));
    psWOptions->panSrcBands[0] = 1;
    psWOptions->panDstBands[0] = 1;
    psWOptions->pfnTransformer = GDALGenImgProjTransform;
    psWOptions->pTransformerArg =
        GDALCreateGenImgProjTransformer( psWOptions->hSrcDS, GDALGetProjectionRef(hSrcDS),
                                         NULL, GDALGetProjectionRef(hSrcDS),
                                         TRUE, 1.0, 0 );
    GDALSetGenImgProjTransformerDstGeoTransform(
        psWOptions->pTransformerArg, adfGeoTransform );
    GDALDatasetH hVRTDS = GDALCreateWarpedVRT( hSrcDS, 1024, 1024,
                                               adfGeoTransform, psWOptions );
    GDALDestroyWarpOptions(psWOptions);
    GDALSetProjection(hVRTDS, GDALGetProjectionRef(hSrcDS));
    GDALSetDescription(hVRTDS, "out.vrt");
    GDALClose(hVRTDS);
    GDALClose(hSrcDS);
    return 0;
}

Running this executable will generate

<VRTDataset rasterXSize="1024" rasterYSize="1024" subClass="VRTWarpedDataset">
  <SRS>PROJCS["NAD27 / UTM zone 11N",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke 1866",6378206.4,294.9786982138982,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG","6267"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26711"]]</SRS>
  <GeoTransform>  4.4072000000000000e+05,  1.1718750000000000e+00,  0.0000000000000000e+00,  3.7513200000000000e+06,  0.0000000000000000e+00, -1.1718750000000000e+00</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTWarpedRasterBand">
    <ColorInterp>Gray</ColorInterp>
  </VRTRasterBand>
  <BlockXSize>512</BlockXSize>
  <BlockYSize>128</BlockYSize>
  <GDALWarpOptions>
    <WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
    <ResampleAlg>NearestNeighbour</ResampleAlg>
    <WorkingDataType>Byte</WorkingDataType>
    <Option name="INIT_DEST">0</Option>
    <SourceDataset relativeToVRT="1">test.tif</SourceDataset>
    <Transformer>
      <GenImgProjTransformer>
        <SrcGeoTransform>440720,0.5859375,0,3751320,0,-0.5859375</SrcGeoTransform>
        <SrcInvGeoTransform>-752162.133333333302,1.70666666666666678,0,6402252.79999999981,0,-1.70666666666666678</SrcInvGeoTransform>
        <DstGeoTransform>440720,1.171875,0,3751320,0,-1.171875</DstGeoTransform>
        <DstInvGeoTransform>-376081.066666666651,0.853333333333333388,0,3201126.39999999991,0,-0.853333333333333388</DstInvGeoTransform>
      </GenImgProjTransformer>
    </Transformer>
    <BandList>
      <BandMapping src="1" dst="1" />
    </BandList>
  </GDALWarpOptions>
</VRTDataset>

One possibility would be to use the return of the GDALCreateOverviewDataset() function as the hSrcDS argument of GDALCreateWarpedVRT(), but GDALCreateOverviewDataset() is C++/internal only.
Another solution, that will work, would be to use the GDALWarp() C function of gdal_utils.h (added in GDAL 2.1), that is the librarified version of gdalwarp itself (if you use it with in-memory datasets, you likely need GDAL 2.2 since I remember having done fixes related to that). And then it will use the overviews, since the source overview selection logic is in gdalwarp_lib.cpp

gdalwarp test.tif test_warped.vrt -of VRT -ts 1024 1024 -overwrite

generates

<VRTDataset rasterXSize="1024" rasterYSize="1024" subClass="VRTWarpedDataset">
  <SRS>PROJCS["NAD27 / UTM zone 11N",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke 1866",6378206.4,294.9786982138982,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG","6267"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26711"]]</SRS>
  <GeoTransform>  4.4072000000000000e+05,  1.1718750000000000e+00,  0.0000000000000000e+00,  3.7513200000000000e+06,  0.0000000000000000e+00, -1.1718750000000000e+00</GeoTransform>
  <Metadata>
    <MDI key="AREA_OR_POINT">Area</MDI>
  </Metadata>
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTWarpedRasterBand">
    <ColorInterp>Gray</ColorInterp>
  </VRTRasterBand>
  <BlockXSize>512</BlockXSize>
  <BlockYSize>128</BlockYSize>
  <GDALWarpOptions>
    <WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
    <ResampleAlg>NearestNeighbour</ResampleAlg>
    <WorkingDataType>Byte</WorkingDataType>
    <Option name="INIT_DEST">0</Option>
    <SourceDataset relativeToVRT="1">test.tif</SourceDataset>
    <OpenOptions>
      <OOI key="OVERVIEW_LEVEL">0</OOI>
    </OpenOptions>
    <Transformer>
      <ApproxTransformer>
        <MaxError>0.125</MaxError>
        <BaseTransformer>
          <GenImgProjTransformer>
            <SrcGeoTransform>440720,1.171875,0,3751320,0,-1.171875</SrcGeoTransform>
            <SrcInvGeoTransform>-376081.066666666651,0.853333333333333388,0,3201126.39999999991,0,-0.853333333333333388</SrcInvGeoTransform>
            <DstGeoTransform>440720,1.171875,0,3751320,0,-1.171875</DstGeoTransform>
            <DstInvGeoTransform>-376081.066666666651,0.853333333333333388,0,3201126.39999999991,0,-0.853333333333333388</DstInvGeoTransform>
          </GenImgProjTransformer>
        </BaseTransformer>
      </ApproxTransformer>
    </Transformer>
    <BandList>
      <BandMapping src="1" dst="1" />
    </BandList>
  </GDALWarpOptions>
</VRTDataset>

Note the addition of

    <OpenOptions>
      <OOI key="OVERVIEW_LEVEL">0</OOI>
    </OpenOptions>

Regarding 'm using a warped VRT because GDAL doesn't surface a non-warped VRT constructor in its C API, that's not actually true. See https://github.com/OSGeo/gdal/blob/trunk/gdal/frmts/vrt/gdal_vrt.h#L84 , but that's a bit involved. You need to create the VRT dataset, add bands, and attach sources to VRT bands. Source overviews will properly be used.
With GDAL 2.1, you can accomplish the same thing more easily with GDALTranslate().
Regular source-based VRT will be faster than warp VRT if you don't need reprojection

sgillies · 2017-10-06T03:34:52Z

I've uploaded wheels made from commit 82beae8 to S3 for testing.

$ aws s3 ls s3://mapbox/rasterio/dist/82beae8/
2017-10-05 19:09:42   25963543 rasterio-1.0a9-cp27-cp27m-macosx_10_6_intel.whl
2017-10-05 19:11:18   52730216 rasterio-1.0a9-cp27-cp27m-manylinux1_x86_64.whl
2017-10-05 19:15:28   52730128 rasterio-1.0a9-cp27-cp27mu-manylinux1_x86_64.whl
2017-10-05 19:19:04   25953080 rasterio-1.0a9-cp34-cp34m-macosx_10_6_intel.whl
2017-10-05 19:20:38   53116657 rasterio-1.0a9-cp34-cp34m-manylinux1_x86_64.whl
2017-10-05 19:26:09   25867888 rasterio-1.0a9-cp35-cp35m-macosx_10_6_intel.whl
2017-10-05 19:27:52   53035929 rasterio-1.0a9-cp35-cp35m-manylinux1_x86_64.whl
2017-10-05 19:32:12   25870078 rasterio-1.0a9-cp36-cp36m-macosx_10_6_intel.whl
2017-10-05 19:33:43   53137374 rasterio-1.0a9-cp36-cp36m-manylinux1_x86_64.whl
2017-10-05 19:36:57    1480895 rasterio-1.0a9.tar.gz

To install the wheel for Python 3.6 on OS X (for example) into a new virtualenv (because this is a pre-release), execute the following command.

pip install https://s3.amazonaws.com/mapbox/rasterio/dist/82beae8/rasterio-1.0a9-cp36-cp36m-macosx_10_6_intel.whl

vincentsarago · 2017-10-06T04:48:27Z

@sgillies Excellent, the wheels work well but the size might be a problem.
The wheel by itself is 53Mb. and when unzipped, the gdal.so is > 100Mb

rouault · 2017-10-06T09:08:42Z

@sgillies Seeing 82beae8 , it looks like you would have managed to make source overviews to be used despite what I thought ?

vincentsarago · 2017-10-06T11:43:01Z

@rouault I can confirm that @sgillies method is working, and fetches overview.

rouault · 2017-10-06T12:45:21Z

Oh ok, I've actually tried this branch and now understand how you use GDAL. Sean's comment "Notice that we read from the VRT at its full resolution" made me take a wrong assumption on how you built the warped VRT. From what I see, you actually build a VRT at its full resolution, but you use GDADatasetRasterIOEx() to read it at lower resolution, and then https://github.com/OSGeo/gdal/blob/8372431040f2b17bb5f6b643c2bfdfcec8d6f630/gdal/frmts/vrt/vrtsourcedrasterband.cpp#L169 kicks in to select the appropriate source overview. Mystery solved

vincentsarago · 2017-10-06T13:18:45Z

Excellent! Thanks @rouault for having a 👀 at it ;-)

sgillies · 2017-10-09T20:19:53Z

Yes, thank you, @rouault! I've used the warped VRT because the non-warped VRT class doesn't expose floating point windows. I'll try to remember to search for or file a GDAL ticket about that.

Sean Gillies added 2 commits October 2, 2017 15:29

Switch to warped vrt internally for boundless reads

e23767b

Update change log, version, fix Window ctor calls

566c6ec

sgillies requested a review from brendan-ward October 3, 2017 18:25

sgillies added this to the 1.0a10 milestone Oct 3, 2017

sgillies added the bug label Oct 3, 2017

Make our temp VRTs larger so we hit the source overviews properly

82beae8

sgillies mentioned this pull request Oct 6, 2017

openjpeg + hdf/netcdf sgillies/frs-wheel-builds#18

Merged

sgillies merged commit 146e2e1 into master Oct 7, 2017

sgillies deleted the boundless-vrt branch October 9, 2017 20:17

markerikson mentioned this pull request Mar 13, 2018

Which GDAL version for compiling geo-data/cesium-terrain-builder#15

Open

vincentsarago mentioned this pull request Dec 5, 2019

Dataset overviews as Rasterio datasets #1504

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Switch to warped vrt internally for boundless reads #1161

Switch to warped vrt internally for boundless reads #1161

sgillies commented Oct 2, 2017 •

edited

Loading

sgillies commented Oct 3, 2017

rouault commented Oct 3, 2017 •

edited

Loading

sgillies commented Oct 6, 2017 •

edited

Loading

vincentsarago commented Oct 6, 2017

rouault commented Oct 6, 2017

vincentsarago commented Oct 6, 2017

rouault commented Oct 6, 2017

vincentsarago commented Oct 6, 2017

sgillies commented Oct 9, 2017

Switch to warped vrt internally for boundless reads #1161

Switch to warped vrt internally for boundless reads #1161

Conversation

sgillies commented Oct 2, 2017 • edited Loading

sgillies commented Oct 3, 2017

rouault commented Oct 3, 2017 • edited Loading

sgillies commented Oct 6, 2017 • edited Loading

vincentsarago commented Oct 6, 2017

rouault commented Oct 6, 2017

vincentsarago commented Oct 6, 2017

rouault commented Oct 6, 2017

vincentsarago commented Oct 6, 2017

sgillies commented Oct 9, 2017

sgillies commented Oct 2, 2017 •

edited

Loading

rouault commented Oct 3, 2017 •

edited

Loading

sgillies commented Oct 6, 2017 •

edited

Loading