Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

VRT support #64

Closed
BoehnkeC opened this issue Jul 28, 2020 · 3 comments
Closed

VRT support #64

BoehnkeC opened this issue Jul 28, 2020 · 3 comments
Labels
enhancement New feature or request

Comments

@BoehnkeC
Copy link

BoehnkeC commented Jul 28, 2020

GDAL offers building of virtual raster datasets, so does rasterio and this gives us the possibility to have it implemented with an interface easy to use.

Consider the following:

import os
import rasterio.shutil

from ukis_pysat import raster as pysat_raster

vrt_dir = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'tests', 'testfiles', 'vrt')

for sample in os.listdir(vrt_dir):
    src = pysat_raster.Image(os.path.join(vrt_dir, sample))
    with rasterio.vrt.WarpedVRT(src.dataset, crs=src.dataset.crs, resampling=rasterio.enums.Resampling.nearest) as vrt:
        rasterio.shutil.copy(vrt, os.path.join(vrt_dir, 'test.vrt'), driver='VRT')

vrt = pysat_raster.Image(os.path.join(vrt_dir, 'test.vrt'))
@BoehnkeC BoehnkeC added enhancement New feature or request good first issue Good for newcomers and removed good first issue Good for newcomers labels Jul 28, 2020
@BoehnkeC
Copy link
Author

BoehnkeC commented Jul 28, 2020

This topic is a bit more complicated. A first draft is provided with this commit.

Having an intial source raster and a bunch of auxiliary rasters, overlapping the source raster, we want to create a mosaic from the auxiliary rasters. This can be clipped to the extent of the source raster later on.

src = Image(source_raster_file)
array, transform = src.mosaic_from_vrt([Image(item) for item in list_of_aux_raster_files])

We call a new method from ukis_pysat.raster, namely mosaic_from_vrt. It is referenced to the source raster object and accepts a list of auxiliary raster objects.

Within mosaic_from_vrt we perform the following steps:

with rasterio.vrt.WarpedVRT(aux_raster_object.dataset, **vrt_options) as vrt:
    target = aux_raster_object.__update_dataset(crs=vrt.crs, transform=vrt.transform)

This produces a VRT dataset, which is basically an XML styled file with transformation information and reference to its original raster object. The VRT dataset can be appended to a list of VRT datasets.

rasterio.merge.merge() accepts a list of datasets and merges them to a resulting array with transformation information.

A few thoughts

  • rasterio.vrt.WarpedVRT keeps a reference of the original dataset in memory: <SourceFilename relativeToVRT="0">/vsimem/5751c324-e285-49ea-ac5a-32b8be93a3fe.</SourceFilename>. Maybe we find a way of referencing files on disk as they are usually located there.
  • Also this comment mentions performance issues of VRT compared to directly warping rasters to a given extent. We have to check that.

@BoehnkeC
Copy link
Author

I compared exploiting the VRT method with the direct use of rasterio.merge.merge, i.e. merging from VRT datasets or building the mosaic from the auxiliary rasters directly. The test revealed that using VRT support is as twice as much slower compared to not using VRT support and building the mosaic right away. Also the memory consumption is a bit higher when using VRT support. The changes were provided with this commit.

This was referenced Jul 30, 2020
@BoehnkeC
Copy link
Author

As they is no need for VRT implementation we decided to drop this issue. Further development will be conducted in #65.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

No branches or pull requests

1 participant