archive-gh-me/rasterio
1
0
Fork 0
mirror of https://github.com/rasterio/rasterio.git synced 2026-02-01 14:34:43 +00:00
Code Issues Projects Releases Wiki Activity

Merge 1.4.3 (#3270)

Browse Source 
* Eliminate boundless reads in merge, new compositing implementation (#3234)

* Eliminate boundless reads in merge, new compositing implementation

Resolves #3228

* Fix nits

* Compare non-zero mean of arrays

* Implement and use Window.align()

Combines the effects of Window.round_offsets and
Window.round_lengths, producing the same effect as gdal_merge.py

* Remove unused math import

* Add more docs about align()

Also increase the fraction used in round_offsets() to be
consistent with align()

* Move align() to a private func, use two existing methods in tests

* Add a test for merging WarpedVRTs (#3237)

Resolves #3196

* Backport of #3217 (#3243)

* Backport of #3217

* Update change log

* Increment GDAL and Python versions for CI (#3244)

* Rewrite _matches() to better support to_authority() and to_epsg() (#3255)

* Rewrite _matches() to better support to_authority() and to_epsg()

Resolves #3239

* Remove list() call and update change log

* Use to_epsg() in is_epsg_code() (#3258)

* Use to_epsg() in is_epsg_code()

Resolves #3248

* Update change log

* Allow MRF compression to surface in properties (#3259)

* Allow MRF compression to surface in properties

Resolves #3256

* Update change log

* Register drivers at most once per process (#3260)

* Register drivers at most once per process

Resolves #3250

* Appease flake8

* Add a note about GDAL_SKIP in the change log

* Support all GDALFillNodata() options in rasterio.fill (#3265)

* Support all GDALFillNodata() options

Resolve #3175.

* Cast values to str and update docs

* Update change log

* Prevent rasterio from trying to open a dataset object (#3266)

Resolves #3105.

* Fix typos discovered by codespell (#3264) (#3267)

* Fix typos discovered by codespell

* crasher



---------

Co-authored-by: Christian Clauss <cclauss@me.com>

* Fix erroneous masking of 0-valued data (#3268)

* Fix erroneous masking of 0-valued data

Resolves #3245

* Add an assertion about data values and update change log

* This is 1.4.3

---------

Co-authored-by: Christian Clauss <cclauss@me.com>
This commit is contained in:
Sean Gillies 2024-12-02 10:43:14 -07:00 committed by GitHub
parent 7a741ee370
commit 24d79d54af
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
19 changed files with 366 additions and 191 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
CHANGES.txt
View File

@ -8,6 +8,35 @@ New features:
- The CRS class has a new, lazily computed, geodetic_crs property (#3218)
1.4.3 (2024-12-02)
------------------
Bug fixes:
- Erroneous masking of 0-valued raster data by boundless, masked reads has been
fixed (#3268).
- If passed a dataset object, rasterio.open() now raises TypeError instead of
proceeding and crashing (#3266).
- All options of GDALFillNodata() are now supported by the method in
rasterio.fill (#3265).
- The flag for GDAL driver registration has been changed to an _env module
attribute. Drivers should only be registered once per process at most
(#3260). A side effect of this is that the GDAL_SKIP configuration option,
which affects format driver registration, only has an effect the first time
a dataset is opened.
- Allow a dataset's compression metadata to surface in profile and compression
properties even if the value isn't present in the Compression enum (#3259).
- A bug that causes CRS.from_wkt().is_epsg_code() to erroneously return False
when an EPSG code is embedded in WKT has been fixed (#3258).
- IAU 2015 has been added to the list of known CRS authorities (#3243).
- A major performance regression in Rasterio's merge tool has been corrected
(#3234).
Other changes:
- CRS._matches() has been rewritten to better support CRS.to_authority() and
CRS.to_epsg() (#3255).
1.4.2 (2024-10-30)
------------------

3
rasterio/__init__.py
View File

@ -22,6 +22,7 @@ if platform.system() == "Windows":
if p and glob.glob(os.path.join(p, "gdal*.dll")):
os.add_dll_directory(os.path.abspath(p))
from rasterio._base import DatasetBase
from rasterio._io import Statistics
from rasterio._vsiopener import _opener_registration
from rasterio._show_versions import show_versions
@ -247,7 +248,7 @@ def open(
or hasattr(fp, "read")
or hasattr(fp, "write")
or isinstance(fp, (os.PathLike, MemoryFile, FilePath))
):
) or isinstance(fp, DatasetBase):
raise TypeError(f"invalid path or file: {fp!r}")
if not isinstance(mode, str):
raise TypeError(f"invalid mode: {mode!r}")

9
rasterio/_base.pyx
View File

@ -981,9 +981,10 @@ cdef class DatasetBase:
def compression(self):
val = self.tags(ns='IMAGE_STRUCTURE').get('COMPRESSION')
if val:
# 'YCbCr JPEG' will be normalized to 'JPEG'
val = val.split(' ')[-1]
return Compression(val)
try:
return Compression(val)
except ValueError:
return val
else:
return None
@ -1037,7 +1038,7 @@ cdef class DatasetBase:
m.update(tiled=self.block_shapes and all(self.width != w for _, w in self.block_shapes))
if self.compression:
m['compress'] = self.compression.name
m['compress'] = getattr(self.compression, "name", self.compression)
if self.interleaving:
m['interleave'] = self.interleaving.name
if self.photometric:

2
rasterio/_env.pxd
View File

@ -6,4 +6,4 @@ cdef class ConfigEnv:
cdef class GDALEnv(ConfigEnv):
cdef public object _have_registered_drivers
pass

111
rasterio/_env.pyx
View File

@ -356,12 +356,14 @@ class PROJDataFinder:
return datadir if os.path.exists(datadir) else None
_have_registered_drivers = False
cdef class GDALEnv(ConfigEnv):
"""Configuration and driver management"""
def __init__(self, **options):
super().__init__(**options)
self._have_registered_drivers = False
def start(self):
CPLPushErrorHandler(<CPLErrorHandler>logging_error_handler)
@ -369,65 +371,16 @@ cdef class GDALEnv(ConfigEnv):
# The outer if statement prevents each thread from acquiring a
# lock when the environment starts, and the inner avoids a
# potential race condition.
if not self._have_registered_drivers:
if not _have_registered_drivers:
with threading.Lock():
if not self._have_registered_drivers:
global _have_registered_drivers
if not _have_registered_drivers:
GDALAllRegister()
OGRRegisterAll()
install_filepath_plugin(filepath_plugin)
if 'GDAL_DATA' in os.environ:
log.debug("GDAL_DATA found in environment.")
self.update_config_options(GDAL_DATA=os.environ['GDAL_DATA'])
else:
path = GDALDataFinder().search_wheel()
if path:
log.debug("GDAL data found in package: path=%r.", path)
self.update_config_options(GDAL_DATA=path)
# See https://github.com/rasterio/rasterio/issues/1631.
elif GDALDataFinder().find_file("gdalvrt.xsd"):
log.debug("GDAL data files are available at built-in paths.")
else:
path = GDALDataFinder().search()
if path:
log.debug("GDAL data found in other locations: path=%r.", path)
self.update_config_options(GDAL_DATA=path)
if 'PROJ_DATA' in os.environ:
# PROJ 9.1+
log.debug("PROJ_DATA found in environment.")
path = os.environ["PROJ_DATA"]
set_proj_data_search_path(path)
elif 'PROJ_LIB' in os.environ:
# PROJ < 9.1
log.debug("PROJ_LIB found in environment.")
path = os.environ["PROJ_LIB"]
set_proj_data_search_path(path)
else:
path = PROJDataFinder().search_wheel()
if path:
log.debug("PROJ data found in package: path=%r.", path)
set_proj_data_search_path(path)
elif PROJDataFinder().has_data():
log.debug("PROJ data files are available at built-in paths.")
else:
path = PROJDataFinder().search()
if path:
log.debug("PROJ data found in other locations: path=%r.", path)
set_proj_data_search_path(path)
if driver_count() == 0:
CPLPopErrorHandler()
raise ValueError("Drivers not registered.")
@ -436,7 +389,57 @@ cdef class GDALEnv(ConfigEnv):
# will acquire a threadlock every time a new environment
# is started rather than just whenever the first thread
# actually makes it this far.
self._have_registered_drivers = True
_have_registered_drivers = True
if 'GDAL_DATA' in os.environ:
log.debug("GDAL_DATA found in environment.")
self.update_config_options(GDAL_DATA=os.environ['GDAL_DATA'])
else:
path = GDALDataFinder().search_wheel()
if path:
log.debug("GDAL data found in package: path=%r.", path)
self.update_config_options(GDAL_DATA=path)
# See https://github.com/rasterio/rasterio/issues/1631.
elif GDALDataFinder().find_file("gdalvrt.xsd"):
log.debug("GDAL data files are available at built-in paths.")
else:
path = GDALDataFinder().search()
if path:
log.debug("GDAL data found in other locations: path=%r.", path)
self.update_config_options(GDAL_DATA=path)
if 'PROJ_DATA' in os.environ:
# PROJ 9.1+
log.debug("PROJ_DATA found in environment.")
path = os.environ["PROJ_DATA"]
set_proj_data_search_path(path)
elif 'PROJ_LIB' in os.environ:
# PROJ < 9.1
log.debug("PROJ_LIB found in environment.")
path = os.environ["PROJ_LIB"]
set_proj_data_search_path(path)
else:
path = PROJDataFinder().search_wheel()
if path:
log.debug("PROJ data found in package: path=%r.", path)
set_proj_data_search_path(path)
elif PROJDataFinder().has_data():
log.debug("PROJ data files are available at built-in paths.")
else:
path = PROJDataFinder().search()
if path:
log.debug("PROJ data found in other locations: path=%r.", path)
set_proj_data_search_path(path)
log.debug("Started GDALEnv: self=%r.", self)

33
rasterio/_fill.pyx
View File

@ -1,4 +1,5 @@
# distutils: language = c++
# cython: c_string_type=unicode, c_string_encoding=utf8
"""Raster fill."""
@ -9,8 +10,13 @@ from rasterio._err cimport exc_wrap_int
from rasterio._io cimport MemoryDataset
def _fillnodata(image, mask, double max_search_distance=100.0,
int smoothing_iterations=0):
def _fillnodata(
image,
mask,
double max_search_distance=100.0,
int smoothing_iterations=0,
**filloptions
):
cdef GDALRasterBandH image_band = NULL
cdef GDALRasterBandH mask_band = NULL
cdef char **alg_options = NULL
@ -27,11 +33,28 @@ def _fillnodata(image, mask, double max_search_distance=100.0,
mask_dataset = MemoryDataset(mask_cast)
mask_band = mask_dataset.band(1)
alg_options = CSLSetNameValue(alg_options, "TEMP_FILE_DRIVER", "MEM")
for k, v in filloptions.items():
k = k.upper()
v = str(v)
alg_options = CSLSetNameValue(alg_options, k, v)
if CSLFindName(alg_options, "TEMP_FILE_DRIVER") < 0:
alg_options = CSLSetNameValue(alg_options, "TEMP_FILE_DRIVER", "MEM")
exc_wrap_int(
GDALFillNodata(image_band, mask_band, max_search_distance, 0,
smoothing_iterations, alg_options, NULL, NULL))
GDALFillNodata(
image_band,
mask_band,
max_search_distance,
0,
smoothing_iterations,
alg_options,
NULL,
NULL
)
)
return np.asarray(image_dataset)
finally:
if image_dataset is not None:
image_dataset.close()

2
rasterio/_io.pyx
View File

@ -695,7 +695,7 @@ cdef class DatasetReaderBase(DatasetBase):
log.debug("Boundless read: self.transform=%r, self.window_transform(window)=%r", self.transform, self.window_transform(window))
mask_vrt_doc = _boundless_vrt_doc(
self, nodata=0,
self, nodata=nodataval,
width=max(self.width, window.width) + 1,
height=max(self.height, window.height) + 1,
transform=self.window_transform(window),

55
rasterio/crs.pyx
View File

@ -20,8 +20,10 @@ not call PROJ functions directly, but invokes them via calls to GDAL's
"""
from collections import defaultdict
from itertools import groupby
import json
import logging
from operator import itemgetter
import pickle
import typing
import warnings
@ -48,6 +50,14 @@ _RE_PROJ_PARAM = re.compile(r"""
""", re.X)
def auth_preference(item):
preferred_order = ["EPSG", "OGC", "ESRI", "IAU_2015"]
conf, auth, val = item
if auth in preferred_order:
return preferred_order.index(auth)
else:
return 100
cdef _safe_osr_release(OGRSpatialReferenceH srs):
"""Wrapper to handle OSR release when NULL."""
@ -140,7 +150,7 @@ cdef class CRS:
"""
try:
return bool(self._epsg)
return bool(self.to_epsg())
except CRSError:
return False
@ -456,11 +466,13 @@ cdef class CRS:
return self._epsg
else:
matches = self._matches(confidence_threshold=confidence_threshold)
if "EPSG" in matches:
self._epsg = int(matches["EPSG"][0])
return self._epsg
else:
if not matches:
return None
else:
for conf, auth, val in matches:
if auth == "EPSG":
return int(val)
return None
def to_authority(self, confidence_threshold=70):
"""Convert to the best match authority name and code.
@ -485,20 +497,17 @@ cdef class CRS:
or None
"""
matches = self._matches(confidence_threshold=confidence_threshold)
# Note: before version 1.2.7 this function only paid attention
# to EPSG as an authority, which is why it takes priority over
# others even if they were a better match.
if "EPSG" in matches:
return "EPSG", matches["EPSG"][0]
elif "OGC" in matches:
return "OGC", matches["OGC"][0]
elif "ESRI" in matches:
return "ESRI", matches["ESRI"][0]
elif "IAU_2015" in matches:
return "IAU_2015", matches["IAU_2015"][0]
if self._epsg is not None:
return ("EPSG", str(self._epsg))
else:
return None
matches = self._matches(confidence_threshold=confidence_threshold)
if not matches:
return None
grouped = groupby(matches, key=itemgetter(0))
for k, group in grouped:
conf, auth, val = sorted(group, key=auth_preference)[0]
return auth, str(val)
def _matches(self, confidence_threshold=70):
"""Find matches in authority files.
@ -516,12 +525,10 @@ cdef class CRS:
cdef int *confidences = NULL
cdef int num_matches = 0
cdef int i = 0
results = defaultdict(list)
cdef list results = []
try:
osr = exc_wrap_pointer(OSRClone(self._osr))
matches = OSRFindMatches(osr, NULL, &num_matches, &confidences)
for i in range(num_matches):
@ -535,12 +542,10 @@ cdef class CRS:
log.debug("returned authority name was null")
if c_code != NULL and c_name != NULL and confidence >= confidence_threshold:
log.debug(
"Matched. confidence=%r, c_code=%r, c_name=%r",
confidence, c_code, c_name)
code = c_code.decode('utf-8')
name = c_name.decode('utf-8')
results[name].append(code)
results.append((confidence, name, code))
return results
finally:

6
rasterio/enums.py
View File

@ -162,7 +162,11 @@ class OverviewResampling(IntEnum):
class Compression(Enum):
"""Available compression algorithms."""
"""Available compression algorithms for GeoTIFFs.
Note that compression options for EXR, MRF, etc are not included
in this enum.
"""
jpeg = 'JPEG'
lzw = 'LZW'
packbits = 'PACKBITS'

14
rasterio/fill.py
View File

@ -9,10 +9,8 @@ from rasterio import dtypes
@ensure_env
def fillnodata(
image,
mask=None,
max_search_distance=100.0,
smoothing_iterations=0):
image, mask=None, max_search_distance=100.0, smoothing_iterations=0, **filloptions
):
"""Fill holes in raster data by interpolation
This algorithm will interpolate values for all designated nodata
@ -48,6 +46,11 @@ def fillnodata(
smoothing_iterations : integer, optional
The number of 3x3 smoothing filter passes to run. The default is
0.
filloptions :
Keyword arguments providing finer control over filling. See
https://gdal.org/en/stable/api/gdal_alg.html. Lowercase option
names and numerical values are allowed. For example:
nodata=0 is a valid keyword argument.
Returns
-------
@ -66,5 +69,4 @@ def fillnodata(
max_search_distance = float(max_search_distance)
smoothing_iterations = int(smoothing_iterations)
return _fillnodata(
image, mask, max_search_distance, smoothing_iterations)
return _fillnodata(image, mask, max_search_distance, smoothing_iterations)

136
rasterio/merge.py
View File

@ -11,9 +11,13 @@ import numbers
import numpy as np
import rasterio
from rasterio.coords import disjoint_bounds
from rasterio.enums import Resampling
from rasterio.errors import MergeError, RasterioDeprecationWarning, RasterioError
from rasterio.errors import (
MergeError,
RasterioDeprecationWarning,
RasterioError,
WindowError,
)
from rasterio.io import DatasetWriter
from rasterio import windows
from rasterio.transform import Affine
@ -396,7 +400,27 @@ def merge(
else:
chunks = [dout_window]
logger.debug("Chunks=%r", chunks)
def _intersect_bounds(bounds1, bounds2, transform):
"""Based on gdal_merge.py."""
int_w = max(bounds1[0], bounds2[0])
int_e = min(bounds1[2], bounds2[2])
if int_w >= int_e:
raise ValueError
if transform.e < 0:
# north up
int_s = max(bounds1[1], bounds2[1])
int_n = min(bounds1[3], bounds2[3])
if int_s >= int_n:
raise ValueError
else:
int_s = min(bounds1[1], bounds2[1])
int_n = max(bounds1[3], bounds2[3])
if int_n >= int_s:
raise ValueError
return int_w, int_s, int_e, int_n
for chunk in chunks:
dst_w, dst_s, dst_e, dst_n = windows.bounds(chunk, output_transform)
@ -404,66 +428,80 @@ def merge(
if inrange:
dest.fill(nodataval)
# From gh-2221
chunk_bounds = windows.bounds(chunk, output_transform)
chunk_transform = windows.transform(chunk, output_transform)
def win_align(window):
"""Equivalent to rounding both offsets and lengths.
This method computes offsets, width, and height that are
useful for compositing arrays into larger arrays and
datasets without seams. It is used by Rasterio's merge
tool and is based on the logic in gdal_merge.py.
Returns
-------
Window
"""
row_off = math.floor(window.row_off + 0.1)
col_off = math.floor(window.col_off + 0.1)
height = math.floor(window.height + 0.5)
width = math.floor(window.width + 0.5)
return windows.Window(col_off, row_off, width, height)
for idx, dataset in enumerate(sources):
with dataset_opener(dataset) as src:
# 0. Precondition checks
# - Check that source is within destination bounds
# - Check that CRS is same
if disjoint_bounds((dst_w, dst_s, dst_e, dst_n), src.bounds):
logger.debug(
"Skipping source: src=%r, bounds=%r",
src,
(dst_w, dst_s, dst_e, dst_n),
)
continue
# Intersect source bounds and tile bounds
if first_crs != src.crs:
raise RasterioError(f"CRS mismatch with source: {dataset}")
# 1. Compute the source window
src_window = windows.from_bounds(
dst_w, dst_s, dst_e, dst_n, src.transform
).round(3)
try:
ibounds = _intersect_bounds(
src.bounds, chunk_bounds, chunk_transform
)
sw = windows.from_bounds(*ibounds, src.transform)
cw = windows.from_bounds(*ibounds, chunk_transform)
except (ValueError, WindowError):
logger.info(
"Skipping source: src=%r, bounds=%r", src, src.bounds
)
continue
temp_shape = (src_count, chunk.height, chunk.width)
cw = win_align(cw)
rows, cols = cw.toslices()
region = dest[:, rows, cols]
temp_src = src.read(
out_shape=temp_shape,
window=src_window,
boundless=True,
masked=True,
if cmath.isnan(nodataval):
region_mask = np.isnan(region)
elif not np.issubdtype(region.dtype, np.integer):
region_mask = np.isclose(region, nodataval)
else:
region_mask = region == nodataval
data = src.read(
out_shape=(src_count, cw.height, cw.width),
indexes=indexes,
masked=True,
window=sw,
resampling=resampling,
)
region = dest[:, :, :]
if cmath.isnan(nodataval):
region_mask = np.isnan(region)
elif not np.issubdtype(region.dtype, np.integer):
region_mask = np.isclose(region, nodataval)
else:
region_mask = region == nodataval
# Ensure common shape, resolving issue #2202.
temp = temp_src[:, : region.shape[1], : region.shape[2]]
temp_mask = np.ma.getmask(temp)
copyto(
region,
temp,
region_mask,
temp_mask,
index=idx,
roff=0,
coff=0,
)
copyto(
region,
data,
region_mask,
data.mask,
index=idx,
roff=cw.row_off,
coff=cw.col_off,
)
if dst:
dst_window = windows.from_bounds(
dst_w, dst_s, dst_e, dst_n, output_transform
).round(3)
dst.write(dest, window=dst_window)
dw = windows.from_bounds(*chunk_bounds, output_transform)
dw = win_align(dw)
dst.write(dest, window=dw)
if dst is None:
if masked:

12
rasterio/windows.py
View File

@ -721,9 +721,8 @@ class Window:
Window
"""
operator = lambda x: int(math.floor(x + 0.5))
width = operator(self.width)
height = operator(self.height)
width = math.floor(self.width + 0.5)
height = math.floor(self.height + 0.5)
return Window(self.col_off, self.row_off, width, height)
def round_shape(self, **kwds):
@ -749,9 +748,8 @@ class Window:
Window
"""
operator = lambda x: int(math.floor(x + 0.001))
row_off = operator(self.row_off)
col_off = operator(self.col_off)
row_off = math.floor(self.row_off + 0.1)
col_off = math.floor(self.col_off + 0.1)
return Window(col_off, row_off, self.width, self.height)
def round(self, ndigits=None):
@ -831,4 +829,4 @@ def subdivide(window, height, width):
row_off += height
col_off = window.col_off
return subwindows
return subwindows

11
tests/test_boundless_read.py
View File

@ -189,3 +189,14 @@ def test_boundless_open_options():
with rasterio.open("tests/data/cogeo.tif", overview_level=2) as src:
data2 = src.read(1, boundless=True)
assert not numpy.array_equal(data1, data2)
def test_issue3245(data):
"""Boundless read of a dataset without a mask should have no masking."""
with rasterio.open(data / "RGB.byte.tif", "r+") as dst:
dst.nodata = None
assert not dst.read(masked=True).mask.any()
data = dst.read(boundless=True, masked=True, window=Window(0, 0, dst.width, dst.height))
assert not data.mask.any()
assert (data >= 0).all()

12
tests/test_crs.py
View File

@ -698,3 +698,15 @@ def test_epsg_4326_ogc_crs84(crs):
"""EPSG:4326 not equivalent to OGC:CRS84."""
assert CRS.from_string("OGC:CRS84") != crs
def test_to_authority():
crs = CRS.from_epsg(32618)
assert crs.to_authority() == ("EPSG", "32618")
def test_is_epsg_code():
crs = CRS.from_wkt(
'PROJCS["RGF93 v1 / Lambert-93",GEOGCS["RGF93 v1",DATUM["Reseau_Geodesique_Francais_1993_v1",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6171"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4171"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",46.5],PARAMETER["central_meridian",3],PARAMETER["standard_parallel_1",49],PARAMETER["standard_parallel_2",44],PARAMETER["false_easting",700000],PARAMETER["false_northing",6600000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2154"]]'
)
assert crs.to_epsg() == 2154
assert crs.is_epsg_code

11
tests/test_env.py
View File

@ -14,7 +14,7 @@ from rasterio import _env
from rasterio._env import del_gdal_config, get_gdal_config, set_gdal_config
from rasterio.env import Env, defenv, delenv, getenv, setenv, ensure_env, ensure_env_credentialled
from rasterio.env import GDALVersion, require_gdal_version
from rasterio.errors import EnvError, RasterioIOError, GDALVersionError
from rasterio.errors import EnvError, GDALVersionError
from rasterio.rio.main import main_group
from rasterio.session import AWSSession, DummySession, OSSSession, SwiftSession, AzureSession
@ -441,13 +441,6 @@ def test_open_with_env(gdalenv):
assert dataset.count == 3
def test_skip_gtiff(gdalenv):
"""De-register GTiff driver, verify that it will not be used."""
with rasterio.Env(GDAL_SKIP='GTiff'):
with pytest.raises(RasterioIOError):
rasterio.open('tests/data/RGB.byte.tif')
@credentials
@pytest.mark.network
def test_s3_open_with_env(gdalenv):
@ -620,7 +613,7 @@ def test_have_registered_drivers():
"""Ensure drivers are only registered once, otherwise each thread will
acquire a threadlock whenever an environment is started."""
with rasterio.Env():
assert rasterio.env.local._env._have_registered_drivers
assert rasterio._env._have_registered_drivers
def test_gdal_cachemax():

93
tests/test_merge.py
View File

@ -12,8 +12,9 @@ import rasterio
from rasterio.merge import merge
from rasterio.crs import CRS
from rasterio.errors import MergeError, RasterioError
from rasterio.vrt import WarpedVRT
from rasterio.warp import aligned_target
from rasterio.windows import subdivide
from rasterio import windows
from .conftest import gdal_version
@ -23,14 +24,31 @@ def test_data_complex(tmp_path):
transform = affine.Affine(30.0, 0.0, 215200.0, 0.0, -30.0, 4397500.0)
t2 = transform * transform.translation(0, 3)
with rasterio.open(tmp_path.joinpath("r2.tif"), 'w', nodata=0, dtype=numpy.complex64, height=2, width=2, count=1,
crs="EPSG:32611", transform=transform) as src:
with rasterio.open(
tmp_path.joinpath("r2.tif"),
"w",
nodata=0,
dtype=numpy.complex64,
height=2,
width=2,
count=1,
crs="EPSG:32611",
transform=transform,
) as src:
src.write(numpy.ones((1, 2, 2)))
with rasterio.open(tmp_path.joinpath("r1.tif"), 'w', nodata=0, dtype=numpy.complex64, height=2, width=2, count=1,
crs="EPSG:32611", transform=t2) as src:
src.write(numpy.ones((1, 2, 2))*2-1j)
with rasterio.open(
tmp_path.joinpath("r1.tif"),
"w",
nodata=0,
dtype=numpy.complex64,
height=2,
width=2,
count=1,
crs="EPSG:32611",
transform=t2,
) as src:
src.write(numpy.ones((1, 2, 2)) * 2 - 1j)
return tmp_path
@ -76,6 +94,7 @@ def test_different_crs(test_data_dir_overlapping):
kwds['crs'] = CRS.from_epsg(3499)
with rasterio.open(test_data_dir_overlapping.joinpath("new.tif"), 'w', **kwds) as ds_out:
ds_out.write(ds_src.read())
with pytest.raises(RasterioError):
result = merge(list(test_data_dir_overlapping.iterdir()))
@ -150,7 +169,7 @@ def test_issue2202(dx, dy):
"/vsis3/copernicus-dem-30m/Copernicus_DSM_COG_10_N48_00_E011_00_DEM/Copernicus_DSM_COG_10_N48_00_E011_00_DEM.tif",
]
]
aux_array, aux_transform = rasterio.merge.merge(datasets=ds, bounds=aoi.bounds)
aux_array, aux_transform = rasterio.merge.merge(ds, bounds=aoi.bounds)
from rasterio.plot import show
show(aux_array)
@ -162,14 +181,16 @@ def test_merge_destination_1(tmp_path):
profile = src.profile
data = src.read()
from rasterio import windows
with rasterio.open(tmp_path.joinpath("test.tif"), "w", **profile) as dst:
for chunk in subdivide(windows.Window(0, 0, dst.width, dst.height), 256, 256):
for chunk in windows.subdivide(
windows.Window(0, 0, dst.width, dst.height), 256, 256
):
chunk_bounds = windows.bounds(chunk, dst.transform)
chunk_arr, chunk_transform = merge([src], bounds=chunk_bounds)
dst_window = windows.from_bounds(*chunk_bounds, dst.transform).round(3)
dst.write(chunk_arr, window=dst_window)
dst_window = windows.from_bounds(*chunk_bounds, dst.transform)
dw = windows.from_bounds(*chunk_bounds, dst.transform)
dw = dw.round_offsets().round_lengths()
dst.write(chunk_arr, window=dw)
with rasterio.open(tmp_path.joinpath("test.tif")) as dst:
result = dst.read()
@ -181,25 +202,31 @@ def test_merge_destination_2(tmp_path):
with rasterio.open("tests/data/RGB.byte.tif") as src:
profile = src.profile
dst_transform, dst_width, dst_height = aligned_target(
src.transform, src.width, src.height, src.res
src.transform,
src.width,
src.height,
src.res,
)
profile.update(transform=dst_transform, width=dst_width, height=dst_height)
data = src.read()
from rasterio import windows
with rasterio.open(tmp_path.joinpath("test.tif"), "w", **profile) as dst:
for chunk in subdivide(windows.Window(0, 0, dst.width, dst.height), 256, 256):
for chunk in windows.subdivide(
windows.Window(0, 0, dst.width, dst.height), 256, 256
):
chunk_bounds = windows.bounds(chunk, dst.transform)
chunk_arr, chunk_transform = merge([src], bounds=chunk_bounds)
dst_window = windows.from_bounds(*chunk_bounds, dst.transform).round(3)
dst.write(chunk_arr, window=dst_window)
dw = windows.from_bounds(*chunk_bounds, dst.transform)
dw = dw.round_offsets().round_lengths()
dst.write(chunk_arr, window=dw)
with rasterio.open(tmp_path.joinpath("test.tif")) as dst:
result = dst.read()
assert result.shape == (3, 719, 792)
assert numpy.allclose(data.mean(), result[:, :-1, :-1].mean())
assert numpy.allclose(
data[data != 0].mean(),
result[result != 0].mean(),
)
@pytest.mark.xfail(gdal_version.at_least("3.8"), reason="Unsolved mask read bug #3070.")
@ -245,3 +272,27 @@ def test_failure_source_transforms(data, matrix):
src.transform = matrix * src.transform
with pytest.raises(MergeError):
merge([src])
def test_merge_warpedvrt(tmp_path):
"""Merge a WarpedVRT into an opened dataset."""
with rasterio.open("tests/data/RGB.byte.tif") as src:
with WarpedVRT(src, crs="EPSG:3857") as vrt:
profile = vrt.profile
data = vrt.read()
profile["driver"] = "GTiff"
with rasterio.open(tmp_path.joinpath("test.tif"), "w", **profile) as dst:
for chunk in windows.subdivide(
windows.Window(0, 0, dst.width, dst.height), 256, 256
):
chunk_bounds = windows.bounds(chunk, dst.transform)
chunk_arr, chunk_transform = merge([vrt], bounds=chunk_bounds)
dst_window = windows.from_bounds(*chunk_bounds, dst.transform)
dw = windows.from_bounds(*chunk_bounds, dst.transform)
dw = dw.round_offsets().round_lengths()
dst.write(chunk_arr, window=dw)
with rasterio.open(tmp_path.joinpath("test.tif")) as dst:
result = dst.read()
assert numpy.allclose(data.mean(), result.mean(), rtol=1e-4)

6
tests/test_open.py
View File

@ -8,6 +8,12 @@ def test_open_bad_path():
rasterio.open(3.14)
def test_open_bad_path_2(path_rgb_byte_tif):
with rasterio.open(path_rgb_byte_tif) as dst:
with pytest.raises(TypeError):
rasterio.open(dst)
def test_open_bad_mode_1():
with pytest.raises(TypeError):
rasterio.open("tests/data/RGB.byte.tif", mode=3.14)

6
tests/test_rio_merge.py
View File

@ -2,7 +2,6 @@
from io import StringIO
import os
import sys
import textwrap
from pathlib import Path
@ -449,8 +448,6 @@ def test_merge_tiny_output_opt(tiffs, runner):
assert data[0][3][0] == 40
@pytest.mark.xfail(sys.version_info < (3,),
reason="Test is sensitive to rounding behavior")
def test_merge_tiny_res_bounds(tiffs, runner):
outputname = str(tiffs.join('merged.tif'))
inputs = [str(x) for x in tiffs.listdir()]
@ -736,7 +733,8 @@ def test_merge_no_gap(tiffs, runner):
with rasterio.open(outputname) as src:
data = src.read(1)
assert data[184, 61] != 0
assert data[183, 61] != 0
assert data[184, 60] != 0
@pytest.mark.parametrize(

6
tests/test_write.py
View File

@ -254,10 +254,10 @@ def test_write_crs_transform_2(tmpdir, monkeypatch):
driver='GTiff', width=100, height=100, count=1,
crs='EPSG:32618',
transform=transform,
dtype=rasterio.ubyte) as s:
s.write(a, indexes=1)
dtype=rasterio.ubyte) as src:
src.write(a, indexes=1)
assert s.crs.to_epsg() == 32618
assert src.crs.to_epsg() == 32618
info = subprocess.check_output(["gdalinfo", name]).decode('utf-8')
assert 'UTM zone 18N' in info
# make sure that pixel size is nearly the same as transform