We've got some new features done and only two outstanding features. Are there any objections to an alpha release tomorrow?

For the wheels, getting this patch https://github.com/OSGeo/gdal/pull/4646 could resolve #2233 .

The Goal

See https://github.com/corteva/rioxarray/issues/246 for details on the origins of this feature.

The goal is to make it as efficient as possible for rasterio to open a Python file-like object. The biggest benefit of a feature like this comes when a library like fsspec is combined with rasterio. fsspec provides a unified interface for accessing files including remote file systems like S3 or GCS and has the ability to cache files locally. It can provide these files as python file-like objects. An example of using fsspec and providing its file to rasterio looks like:

import rasterio
import fsspec

with fsspec.open('https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/13/S/DV/2020/4/S2B_13SDV_20200428_0_L2A/TCI.tif') as fobj:
    with rasterio.open(fobj) as src:
        print(src.profile)

The Current Issue

Currently, rasterio uses a MemoryFile, a GDAL supported feature that reads all the necessary data into memory and accesses the data from there instead of disk. This is great, except that it requires reading all of the data from the file before we need it. This means a lot of time (especially for remote resources) and a lot of memory. With the above example and using the rasterio master branch it takes about ~21s on my laptop to do the above.

The Solution - Attempt 2

The current version of this pull request uses GDAL's plugin interfaces to add all necessary hooks (callbacks) for GDAL to communicate with rasterio and the file-like objects that you provide to it. This comes as two parts: the virtual filesystem and the individual file handling. Right now the filesystem is all kept in a single python dictionary (virtual filename -> file wrapper object). The file handling is fairly basic and doesn't need to be too complicated as the basic FILE* operations map to file-like objects (ex read, seek, tell, etc).

I haven't done anything for writing or for multi-range reads.

TODO

• [x] Switch to GDAL Plugin system. I didn't know about this when I started and it could mean moving all the .cpp stuff into Cython which would be great for maintainability. This will be "Attempt 2".
• [ ] Figure out multi-threading behavior and multi-process behavior. This is mostly related to how the file system treats a file when rasterio is using it. Also, determines if/when we need locking. This includes the GIL. For example, if GDAL is making its own threads (outside of python), then we may need to explicitly create the thread-local GIL.
• [x] Add tests
• [x] Add documentation
• [x] Add docstrings
• [x] Answer below questions

Questions

1. Will this replace all rasterio handling of non-string inputs (objects)? Probably not, but in that case when is MemoryFile used instead? Do users get a kwarg to choose? Answer: It will replace reading of file-like objects.
2. Should this be limited to fsspec objects or any file-like object (handling any missing methods gracefully)? Limiting to fsspec objects may let us assume some things. Answer: All file-like objects.
3. Are users allowed to provide a filename for the open file object to help GDAL know how to refer to it? If so, what happens if a user provides two different objects with the same name? Answer: Yes, but only to match the existing interface of MemoryFile. Conflicts have an undefined behavior.
4. Can we support multi-range reads in a semi-standard way? Answer: No. Ignore it for now.

I've been away from Windows for over a decade. Anyone have experience making Windows binaries with Appveyor as described in https://packaging.python.org/en/latest/appveyor.html?

/cc @brendan-ward

Towards closing #175.

Usage, for now, like this:

$ rio calc "0.10*{1} + 125" tests/data/shade.tif out.tif

Results in a new hillshade with shade values of 125 instead of 0.

Expected behavior and actual behavior.

The Satpy project has a CI environment that installs various dependencies from their "unstable" sources (nightly builds, git repositories, etc). Recently (in the last week probably), I've started getting this exception:

 /usr/share/miniconda3/envs/test-environment/lib/python3.9/site-packages/trollimage/xrimage.py:58: in <module>
    from rasterio.windows import Window
/usr/share/miniconda3/envs/test-environment/lib/python3.9/site-packages/rasterio/__init__.py:11: in <module>
    with rasterio._loading.add_gdal_dll_directories():
E   AttributeError: partially initialized module 'rasterio' has no attribute '_loading' (most likely due to a circular import)

Steps to reproduce the problem.

I'm still trying to reproduce it locally. Theoretically just doing the import above should cause it.

Failing job: https://github.com/pytroll/satpy/actions/runs/3055333856/jobs/4956174627

Environment Information

• Output from: rio --show-versions or python -c "import rasterio; rasterio.show_versions()"

• rasterio version (python -c "import rasterio; print(rasterio.__version__)")

    Created wheel for rasterio: filename=rasterio-1.4.dev0-cp39-cp39-linux_x86_64.whl size=1474741 sha256=ebce139195b79b11d1baa440541745a20a689ef45d6fbcd3af0794544d17e758
  

• GDAL version (python -c "import rasterio; print(rasterio.__gdal_version__)"): 3.5.2

• Python version (python -c "import sys; print(sys.version.replace('\n', ' '))"): 3.9

• Operation System Information (python -c "import platform; print(platform.platform())"): Ubuntu (ubuntu-latest github action)

Installation Method

Conda-forge base packages but pip install various packages with pip from git/nightly builds.

Hello,

When working on Linux (Windows is fine), I have an aliased result on overlapping areas when merging datasets with max and min methods. first method is OK.

The working area:

Here is the correct screenshot (first method):

And here the broken one (both min and max are a like, I didn't test other methods):

When zooming to the pixel, this is what we got:

My code is:

# From the zip uploaded hereafter
paths = ["14_20200228T143729_S2_T18HXC_L1C_175923_CLIP.tif", "14_20200228T143729_S2_T18HXD_L1C_175923_CLIP.tif"]
method = "max"  # min or first

# Open datasets for merging
datasets = []
try:
    for path in paths:
        datasets.append(rasterio.open(str(path)))

    # Merge all datasets
    merged_array, merged_transform = merge.merge(datasets, method=method)
    merged_meta = datasets[0].meta.copy()
    merged_meta.update(
        {
            "driver": "GTiff",
            "height": merged_array.shape[1],
            "width": merged_array.shape[2],
            "transform": merged_transform,
        }
    )
finally:
    # Close all datasets
    for dataset in datasets:
        dataset.close()

# Write on disk
with rasterio.open("test_min2.tif", "w", **merged_meta) as dst:
    dst.write(merged_array)

It does not seem to happen all the time, but here are my datasets: MOSAIC.zip

And my results: results.zip

Operating system

Linux (through docker), not reproducible on windows !

Rasterio version and provenance

Rasterio 1.2.6 via pip (both on windows and linux)

There are a number of places where a rasterio operation takes as input a raster or raster band, OR a numpy ndarray, transform, and sometimes CRS. This complicates the function signatures a bit, as it means having parameters that are meaningful in some cases but not others, and it lacks a bit of cohesion between related entities that may be helpful in many cases.

Perhaps it is time to create a lightweight object (tuple or class) that wraps these together, so that we can pass just one thing around and have it be more equivalent to passing a raster / raster band?

Something like the Raster object in @perrygeo 's python-rasterstats package, but not wrapping rasters / bands; only ndarrays.

This surfaced in #578, among others.

Closes https://github.com/mapbox/rasterio/issues/516 Closes https://github.com/mapbox/rasterio/issues/794

Remove deprecated features and API's. See docs/migrating-to-v1.rst for a general description.

Let's do this before the end of October.

New features include:

• VSI cache clearing #1078
• API for driver metadata (like https://github.com/Toblerity/Fiona/pull/950)
• RPC support #1922
• More: https://github.com/mapbox/rasterio/milestone/95?closed=1

Between 0.29.0 and 0.30.0, functionality was introduced to specifically support the virtual file system functinality vsigzip, vsizip, and vsitar. As can be seen in this code, along with explicit support for some virtual file systems, some functionality for other virtual file system support that was previously and unintentionally exposed was explicitly disallowed, and calls that had taken advantage of that GDAL functionality began throwing a "VFS scheme {0} unknown" exception upon upgrade of rasterio.

Part of the GDAL functionality that became inaccessible through rasterio was accessing rasters over HTTP with vsicurl. That feature in particular was useful to myself and others, since we were relying on it to do windowed reads of large rasters over the network, without pulling down the whole file. In this way we can have multiple processes pulling byte ranges of a GeoTiff off of S3 (which supports the byte range HTTP request header), as a mechanism of doing tiling work using rasterio and PySpark (see the OpenAerialMap tiler code here...this "chunking" step had some problems dealing with very large rasters so if you're interested in using these techniques make sure to look to the master branch's version and any unmerged PR's).

I started a discussion on a merged pull request, which is the place I found the most definitive language about why vsicurl was no longer supported. It was smartly recommended that we start a separate ticket instead of talking on a merged PR, and this ticket can serve this purpose. Please read the prior discussion as if they were the first comments on this issue.

/cc @sgillies @perrygeo @mojodna

Lately, we have encountered a strange bug in Terracotta. It basically always leads to errors like these:

(from https://github.com/DHI-GRAS/terracotta/issues/139)

Traceback (most recent call last):
  File "rasterio/_base.pyx", line 213, in rasterio._base.DatasetBase.__init__
  File "rasterio/_shim.pyx", line 64, in rasterio._shim.open_dataset
  File "rasterio/_err.pyx", line 205, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_OpenFailedError: '/vsis3/italy-composite/rasters/italy_2018_red.tif' not recognized as a supported file format.

or

(from https://github.com/DHI-GRAS/terracotta/issues/10#issuecomment-489826885)

Traceback (most recent call last):
  File "rasterio/_io.pyx", line 698, in rasterio._io.DatasetReaderBase._read
  File "rasterio/shim_rasterioex.pxi", line 133, in rasterio._shim.io_multi_band
  File "rasterio/_err.pyx", line 182, in rasterio._err.exc_wrap_int
rasterio._err.CPLE_AppDefinedError: IReadBlock failed at X offset 0, Y offset 0: '/vsis3/bucket/prefix/tile330.tif' does not exist in the file system, and is not recognized as a supported dataset name.

The errors occur on different versions of rasterio, although anecdotally it wasn't a problem pre-1.0.15. It also seems to occur both during rasterio.open, and when actually reading tiles via WarpedVRT.read.

The problem is that we have only observed it with huge raster files, and we haven't been able to reproduce this reliably, or in a way where I could share it with you.

I am opening this issue to see if you have some intuition why the problem might occur. If you think it just adds noise, please feel free to close it.

For example, gdal_translate uses a VRT under the hood to gather changes and options and then relies on GDALCreateCopy to materialize those changes and produce an output file: https://github.com/OSGeo/gdal/blob/master/apps/gdal_translate_lib.cpp#L2511.

Rasterio's merge tool could similarly take advantage of an internal VRT. Our existing merge tool produces very large numpy arrays as its only possible output. An example of the direction we might take is shown in #2699.

Rasterio stack CLI program could also use a VRT. Stacking by itself wouldn't need to use a lot of resources, it would only be when materializing hefty stacks to a file that the deferred work would be done.

And of course, Rasterio's convert CLI program could use a VRT in the same way gdal_translate does.

There's no immediate action here. I want to show that #2699 isn't entirely out of the blue and is part of a theme I've been thinking about for a while: Rasterio making virtual datasets using Python XML libs and GDAL reading those datasets.

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 #507 and #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 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.

Resolves #2573

(Following https://github.com/rasterio/rasterio/issues/2696)

I have been looking into the implementation of rasterio.merge.merge to "prefetch" common bounds from a group of rasters. Once working on my implementation (since the function isn't currently too modular) I realized that I was getting errors when calculating common bounds of tiffs.

After some investigation, it turns out this was because the CRS of each tif was different, resulting in bounds given in different units. Before implementing my own solution, I looked at the merge implementation again but found no consideration for this. In the end, I implemented a solution based on rasterio.warp.transform_bounds.

I currently implemented this solution with a separate function since it may be a backward breaking change and I'm not sure whether that is something that would be appreciated. I replicated all existing tests for the new merge function (with support for multiple CRSs) but still need to add tests to check the feature specifically (i.e. actively try to merge tiffs with different CRSs).

Discussed in https://github.com/rasterio/rasterio/discussions/2682