Instead of `band_number = np.arange(count) + 1` on line 21, try `band_number = range(1, count + 1)` I didn't test your code specifically, but rasterio won't read a dataset when passed the band indexes

I just tried on Windows and Linux and your code ran fine for me (using the test data you attached). The only other thing I can think of is I'm running rasterio/gdal installed from the conda-forge repo

Sorry Matt, I don't know. Even Rouault wrote a bit about the issue a number of years ago, but the re-writing without re-compressing was left in the "Ideas for later..." section.

Your code runs fine for me Matt. Update your rasterio version? I'm running: $ rio --version 1.2.10 $ gdalinfo --version GDAL 3.4.0, released 2021/11/04 You should note that by re-writing the file, you

Replacing your FTP code with direct access using the /vsicurl/ virtual filesystem: import rasterio vrt = ''' <VRTDataset rasterXSize="6935" rasterYSize="3351"> <SRS>EPSG:4326</SRS> <GeoTransform>-124.

You could use a raw raster VRT combined with the /vsitar and /vsigzip virtual filesystems (you could probably use /vsicurl as well to replace your FTP code, but I couldn't get that to work immediately

You don't need to pre-populate a dataset with a massive empty array, you can just open it in write mode and write chunks of data via windowed writing: profile = { 'driver': 'GTiff', 'tiled': True, 'he

rasterio.transform.from_origin(west, north, xsize, ysize) Return an Affine transformation given upper left and pixel sizes. https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#raster

You can also use json: import json poly = json.loads(obj.job_loc.geojson)

rasterio.mask.mask returns a two element tuple - element 0 is the array, element 1 is the transform. Try something like (untested): # Mask Raster masked_raster, transform = rasterio.mask.mask(file, po

GDALs Sentinel2 driver exposes the data as subdatasets, with (sort of) one for each resolution - https://gdal.org/drivers/raster/sentinel2.html#level-2a To access the subdataset that contains the 10m

Have you tried keeping your windows all the same size (including the edge windows) and specifying `boundless=True` and a `fill_value=some appropriate value` to fill out your windows in your ` `dataset

rasterio.transform.array_bounds should do the trick - https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#rasterio.transform.array_bounds

Try going the other way - multiply the resolution by a scale factor (i.e. 3 in your case) and divide the dimensions - see example gist

No. Exactly as I show in my code example. No need to write it out again, as you are just updating the georeferencing with the dataset opened in update mode ("r+"). The dataset is saved automatically w

To override the projection with rasterio: import rasterio from rasterio.crs import CRS in_tif = 'path/to.tif' with rasterio.open(in_tif, 'r+') as src: src.crs = CRS.from_epsg(32629)

Is your matrix of values a numpy array? You can easily create a GeoTIFF in rasterio if you have an array of data and the basic metadata (georeferencing, dimensions, band count, data type). https://ras

How did you install rasterio? If you installed via Anaconda / conda using the default channel, your gdal may lack bigtiff support. https://github.com/ContinuumIO/anaconda-issues/issues/9887

You could use a context manager to clean it up automatically: from contextlib import contextmanager import rasterio from rasterio import MemoryFile @contextmanager def mem_raster(data, **profile): wit

I've put a couple of methods in my answer to your SO question, by expanding the block windows and also for arbitrary windows (unrelated to internal blocks/tiles).