VRT form of merging #2699

sgillies · 2023-01-02T23:40:47Z

There's a lot of overlap between rasterio's merge and GDAL's GDALBuildVRT utility. I believe having two different ways to do the same kind of thing is a step backwards for this project. As an alternative to #2573, which wraps GDALBuildVRT, I've reimplemented the GDAL utility in Python, building on existing code in _boundless_vrt_doc.

Advantages over using GDALBuildVRT:

It's compatible with Rasterio dataset objects. GDALBuildVRT only understands GDAL filenames.
It returns an XML string, making temporary files on the filesystem unnecessary.
It's Python. Easy to read and modify. GDALBuildVRT is C++. Harder to read and modify.
It shares code with _boundless_vrt_doc. We can continue to refactor and develop a better virtual dataset story throughout Rasterio.
It shares code with merge. The functions could be refactored so that merge uses virtual_merge, which also has the potential to resolve other issues in our tracker (like Improve performance of rio-merge #507 and [WIP] Tiled merge #2221).

The name of this function and the module it lives in are undecided.

This can totally replace gdal.BuildVRT for OSMnx (in #2573) and doesn't require wrapping GDALBuildVRT. Just because a Rasterio user wants something that GDALX does doesn't mean we need to literally expose GDALX from Rasterio. GDAL's API design can be improved upon and Rasterio is the project to do those improvements.

Resolves #2573

sgillies · 2023-01-03T00:16:24Z

@snowman2 the code here is basically correct! The XML produced by virtual_merge(glob("tests/data/rgb?.tif")) becomes

What do you think? I want Rasterio to aim high when it comes to API design. It's pragmatic to wrap gdalbuildvrt, but I think we can do better.

snowman2 · 2023-01-03T14:45:29Z

@sgillies, your logic makes sense to me - thanks for working on this 👍

rasterio/merge.py

gboeing · 2023-01-05T18:52:10Z

This can totally replace gdal.BuildVRT for OSMnx (in #2573) and doesn't require wrapping GDALBuildVRT. Just because a Rasterio user wants something that GDAL.X does doesn't mean we need to literally expose GDAL.X from Rasterio. GDAL's API design can be improved upon and Rasterio is the project to do those improvements.

@sgillies that reasoning makes perfect sense. Thanks!

sgillies · 2023-01-06T17:00:58Z

@gboeing in the meanwhile, you have the option of copying the new virtual_merge function into osmnx if you want. I don't imagine the interface or functionality will change much between now and 1.4.0.

snowman2 · 2023-01-06T17:41:12Z

@sgillies, thoughts about adding more tests? Could start by copying some from #2698.

I am thinking it would be good to verify that the bounds/transform/data of the merged VRT are as expected. Also probably a good idea to test merging multiple files.

If rasters with different CRS are passed into this function, should it raise an exception? Should we test that?

rasterio/merge.py

snowman2 · 2023-01-06T17:43:34Z

rasterio/merge.py

+        Whether to adjust output image bounds so that pixel coordinates
+        are integer multiples of pixel size, matching the ``-tap``
+        options of GDAL utilities.  Default: False.
+    dst_path : str or PathLike, optional


Is dst_path used?

Is dst_kwds used?

These are left over from the original merge. TODO: remove them.

sgillies · 2023-01-06T18:07:12Z

@snowman2 @groutr (adding you since you're a key merge contributor) if we have a virtual merging method like this one, what would you think about it returning a "Virtual Dataset" object instead of XML or an infoset?

It could be very lightweight, like our existing MemoryDataset. MemoryDataset is a thin wrapper around a "/vsimem/..." dataset. A virtual dataset could be a thin wrapper around an XML bytestring. The usage could be a lot like MemoryDataset's. I've got an additional method in mind: a VirtualDataset.materialize("example.tif") that would call GDALCopyCreate to make a non-virtual dataset.

I'll throw a quick example together.

rasterio/merge.py

groutr · 2023-01-06T18:41:12Z

@snowman2 @groutr (adding you since you're a key merge contributor) if we have a virtual merging method like this one, what would you think about it returning a "Virtual Dataset" object instead of XML or an infoset?

It could be very lightweight, like our existing MemoryDataset. MemoryDataset is a thin wrapper around a "/vsimem/..." dataset. A virtual dataset could be a thin wrapper around an XML bytestring. The usage could be a lot like MemoryDataset's. I've got an additional method in mind: a VirtualDataset.materialize("example.tif") that would call GDALCopyCreate to make a non-virtual dataset.

I'll throw a quick example together.

I like the idea, though I would suggest choosing a method name that is easier to type like VirtualDataset.compute("example.tif") (to mirror the dask ecosystem terminology for realizing a result) or VirtualDataset.export("example.tif") or VirtualDataset.output("example.tif"). Another option that might fit well with rasterio is VirtualDataset.to_dataset("example.tif") (which seems like a very natural way to expect turning a virtual dataset into an actual dataset).

snowman2 · 2023-01-06T18:46:47Z

A virtual dataset could be a thin wrapper around an XML bytestring.

I like that idea.

groutr · 2023-01-06T19:39:21Z

A virtual dataset could be a thin wrapper around an XML bytestring.

We might consider adding internal validation against the VRT schema (https://raw.githubusercontent.com/OSGeo/gdal/master/data/gdalvrt.xsd). xmlschema seems like a nice package for validation.

rasterio/_io.pyx

sgillies · 2023-01-10T00:21:34Z

rasterio/vrt.py

+
+        """
+        # To avoid a circular import.
+        from rasterio.io import MemoryFile


Would it make sense to put VirtualDataset with the rest of the Dataset objects in rasterio.io?

rasterio/vrt.py

snowman2 · 2023-01-10T04:13:32Z

rasterio/vrt.py

+
+        kwargs : dict
+            Optional dataset opening options to be passed to
+            rasterio.open().


Suggested change

rasterio.open().

:func:`rasterio.open`.

rasterio/vrt.py

groutr · 2023-01-10T15:49:42Z

rasterio/vrt.py

+        from rasterio.io import MemoryFile
+
+        tree = ET.fromstring(text)
+        doc = ET.tostring(tree)


I'm not understanding why we are parsing an xml string into an elementtree, only to write it back out to an xml string. Can we not use text as is?

@groutr I didn't have any intention other than 1) very basic validation (is it correct XML?) and 2) convert from bytes/unicode to just bytes.

I think we may well want to do schema validation. But for now, this is just a sketch.

snowman2 · 2023-01-10T20:47:52Z

@sgillies, I am thinking that it may make sense to move the VirtualDataset to a separate PR. I think that will help focus the discussion. Thoughts?

groutr · 2023-01-11T05:46:55Z

Is the goal of the VirtualDataset to replicate the functionality of a rasterio Dataset, permit construction and then be able to render a vrt xml or is it simply serving as a VRT wrapper for MemoryFile?

Instead of GDALBuildVRT Resolves #2573

Allow rasterio.open to take XML bytes as well as unicode strings.

snowman2 · 2023-06-08T01:55:30Z

rasterio/merge.py

+    else:
+
+        @contextmanager
+        def nullcontext(obj):


Since this could be confused with contextlib.nullcontext, I propose renaming this to something like:

Suggested change

def nullcontext(obj):

def dataset_opener(obj):

And removing dataset_opener = nullcontext below.

Given that rasterio supports python 3.9+, why not use contextlib.nullcontext directly?

The function defined here behaves differently from contextlib.nullcontext. nullcontext returns None whereas this function returns the file handle.

contextlib.nullcontext will return whatever you pass as the argument (the default is None).
https://docs.python.org/3/library/contextlib.html#contextlib.nullcontext

Nice - then I agree with your suggestion to use contextlib.nullcontext.

@groutr great suggestion! I will make it so.

snowman2 · 2023-06-08T01:57:50Z

rasterio/merge.py

+
+    ET.SubElement(vrtdataset, "SRS").text = crs_wkt if crs_wkt else ""
+    ET.SubElement(vrtdataset, "GeoTransform").text = ",".join(
+        [str(v) for v in output_transform.to_gdal()]


I think this should work as well:

Suggested change

[str(v) for v in output_transform.to_gdal()]

str(v) for v in output_transform.to_gdal()

snowman2 · 2023-06-08T02:00:06Z

rasterio/merge.py

+            nodataval = nodata
+        else:
+            warnings.warn(
+                "The nodata value, %s, is beyond the valid "


Suggested change

"The nodata value, %s, is beyond the valid "

f"The nodata value, {nodataval}, is beyond the valid "

snowman2 · 2023-06-08T02:00:23Z

rasterio/merge.py

+        else:
+            warnings.warn(
+                "The nodata value, %s, is beyond the valid "
+                "range of the chosen data type, %s. Consider overriding it "


Suggested change

"range of the chosen data type, %s. Consider overriding it "

f"range of the chosen data type, {dt}. Consider overriding it "

snowman2 · 2023-06-08T02:00:37Z

rasterio/merge.py

+            warnings.warn(
+                "The nodata value, %s, is beyond the valid "
+                "range of the chosen data type, %s. Consider overriding it "
+                "using the --nodata option for better results." % (nodataval, dt)


Suggested change

"using the --nodata option for better results." % (nodataval, dt)

"using the --nodata option for better results."

sgillies · 2023-06-12T21:06:31Z

Thanks for all the suggestions @groutr @snowman2 ! I'm going to accept them as soon as I'm back from my Wisconsin trail running vacation 😆

sgillies · 2023-09-01T19:42:18Z

Deferring until 1.5.0. There's a lot to this and it's not ready.

sgillies added the vrt label Jan 2, 2023

sgillies self-assigned this Jan 2, 2023

sgillies mentioned this pull request Jan 3, 2023

Using more intermediate VRTs in higher-level tools #2700

Open

This was referenced Jan 3, 2023

ENH: Add rasterio.vrt.BuildVRT #2573

Closed

ENH: Add BuildVRT #2698

Closed

rouault reviewed Jan 4, 2023

View reviewed changes

rasterio/merge.py Outdated Show resolved Hide resolved

sgillies commented Jan 4, 2023

View reviewed changes

rasterio/merge.py Show resolved Hide resolved

sgillies marked this pull request as ready for review January 4, 2023 22:55

sgillies marked this pull request as draft January 4, 2023 22:55

gboeing mentioned this pull request Jan 5, 2023

drop GDAL dependency gboeing/osmnx#916

Open

sgillies marked this pull request as ready for review January 6, 2023 17:01

snowman2 reviewed Jan 6, 2023

View reviewed changes

rasterio/merge.py Show resolved Hide resolved

snowman2 reviewed Jan 6, 2023

View reviewed changes

groutr reviewed Jan 6, 2023

View reviewed changes

rasterio/merge.py Outdated Show resolved Hide resolved

sgillies commented Jan 10, 2023

View reviewed changes

rasterio/_io.pyx Show resolved Hide resolved

sgillies commented Jan 10, 2023

View reviewed changes

rasterio/vrt.py Outdated Show resolved Hide resolved

snowman2 reviewed Jan 10, 2023

View reviewed changes

rasterio/vrt.py Outdated Show resolved Hide resolved

snowman2 reviewed Jan 10, 2023

View reviewed changes

rasterio/vrt.py Outdated Show resolved Hide resolved

groutr reviewed Jan 10, 2023

View reviewed changes

sgillies added this to the 1.4.0 milestone Apr 27, 2023

sgillies added 5 commits May 25, 2023 20:17

VRT form of merging

ec8bb35

Instead of GDALBuildVRT Resolves #2573

Treat XML as UTF-8 encoded bytes

052057b

Allow rasterio.open to take XML bytes as well as unicode strings.

Prototype of VirtualDataset class

2d03100

Ensure MemoryFile is removed

bf1ad79

Back out VirtualDataset

54b21ee

sgillies force-pushed the merge-vrt branch from 1002959 to 54b21ee Compare May 26, 2023 02:22

sgillies added 2 commits May 26, 2023 16:58

Remove unused imports

8c96ff3

Add resampling

14192e3

snowman2 reviewed Jun 8, 2023

View reviewed changes

sgillies removed this from the 1.4.0 milestone Sep 1, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

VRT form of merging #2699

VRT form of merging #2699

sgillies commented Jan 2, 2023 •

edited

sgillies commented Jan 3, 2023

snowman2 commented Jan 3, 2023

gboeing commented Jan 5, 2023

sgillies commented Jan 6, 2023

snowman2 commented Jan 6, 2023

snowman2 Jan 6, 2023

snowman2 Jan 6, 2023

sgillies Jun 12, 2023

sgillies commented Jan 6, 2023

groutr commented Jan 6, 2023

snowman2 commented Jan 6, 2023

groutr commented Jan 6, 2023

sgillies Jan 10, 2023

groutr Jan 10, 2023

snowman2 Jan 10, 2023

groutr Jan 10, 2023 •

edited

sgillies Jan 10, 2023 •

edited

snowman2 commented Jan 10, 2023

groutr commented Jan 11, 2023

snowman2 Jun 8, 2023

groutr Jun 12, 2023

snowman2 Jun 12, 2023

groutr Jun 12, 2023 •

edited

snowman2 Jun 12, 2023

sgillies Jun 12, 2023

snowman2 Jun 8, 2023

snowman2 Jun 8, 2023

snowman2 Jun 8, 2023

snowman2 Jun 8, 2023

sgillies commented Jun 12, 2023

sgillies commented Sep 1, 2023 •

edited

	[str(v) for v in output_transform.to_gdal()]
	str(v) for v in output_transform.to_gdal()

	"The nodata value, %s, is beyond the valid "
	f"The nodata value, {nodataval}, is beyond the valid "

	"range of the chosen data type, %s. Consider overriding it "
	f"range of the chosen data type, {dt}. Consider overriding it "

	"using the --nodata option for better results." % (nodataval, dt)
	"using the --nodata option for better results."

VRT form of merging #2699

Are you sure you want to change the base?

VRT form of merging #2699

Conversation

sgillies commented Jan 2, 2023 • edited

sgillies commented Jan 3, 2023

snowman2 commented Jan 3, 2023

gboeing commented Jan 5, 2023

sgillies commented Jan 6, 2023

snowman2 commented Jan 6, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sgillies commented Jan 6, 2023

groutr commented Jan 6, 2023

snowman2 commented Jan 6, 2023

groutr commented Jan 6, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

groutr Jan 10, 2023 • edited

Choose a reason for hiding this comment

sgillies Jan 10, 2023 • edited

Choose a reason for hiding this comment

snowman2 commented Jan 10, 2023

groutr commented Jan 11, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

groutr Jun 12, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sgillies commented Jun 12, 2023

sgillies commented Sep 1, 2023 • edited

sgillies commented Jan 2, 2023 •

edited

groutr Jan 10, 2023 •

edited

sgillies Jan 10, 2023 •

edited

groutr Jun 12, 2023 •

edited

sgillies commented Sep 1, 2023 •

edited