Here's a very naive way using rasterio.shutil.copy to generate VRT XML for each source raster, combine them and output a stacked VRT (equivalent to gdalbuildvrt -separate).
By naive, I mean it does not check any of the properties and assumes rasters are all the same shape, extent, pixel size, data type, projection etc.
import xml.etree.ElementTree as ET
If you're looking for a QGIS solution, you can check the r.to.vect GRASS function under the QGIS Processing Toolbox.
These links below might help:
Converting raster to vector by generating center lines?
You're writing a floating point array as an unsigned 8 bit integer. Try
src = rasterio.open('raster2.tif')
profile = src.profile
with rio.open('test.tif', 'w', **profile) as dst:
dst.write(new_raster, 1) # No need to cast to int, as profile is already float32
or if that doesn't work, enforce the casting using .astype(rio.float32) instead of your ...
If you're ultimately working toward a polygon (circle) intersection with your raster, there are built-in methods for calculating a numpy mask, given a polygon.
You can use this question to see how Shapely can be used to get a circle of defined radius from a centre point (you will need to convert to a projected coordinate system).
Then, use rasterio.features....
So, I got it working using rasterio. The problem was that the Affine pivot parameter is actually in pixel coordinates, and I was passing real world coordinates. So few lines of code later, this is rotating as exected!
from optparse import OptionParser
from affine import Affine # For easly manipulation of affine matrix
Here's a non geocube way using rasterio and fiona (could also use the higher level geopandas which uses fiona under the hood). This
loops through small chunks (aka a rasterio.windows.Window) and the features in your polygon dataset,
rasterizes the features just in the small chunk area,
sums the chunks,
writes each summed chunk to the appropriate part of ...
# Look at this simple code:
path = r'inputpath_directory'
# for Loop
bands = 
for i in os.listdir(path):
print(bands) # a list of bands
# To save a raster using gdal:
# Look At this code:
import numpy as np
import matplotlib.pyplot as plt
# open the bands:
B4 = rasterio.open(inputpath_name)
B5 = rasterio.open(inputpath_name)
# First calculate your NDVI with this function:
def ndvi(red, nir):
red = B4.read(1)
nir = B5.read(1)
red = red.astype('float64')
nir = nir.astype('float64')
This should work with what you require. Read about rasterio write here
d = ['Blue', 'Red', 'Green']
test = 'test.tif'
folder = Path('Path/to/My/Folder')
with rasterio.open("test.tif") as src:
for f in folder.glob("*.tif"):
if any(color in f.name for color in d):
with rasterio.open(f) as med:
new_raster = ...
I would recommend using @radouxju approach to calculating NDVI in this answer. Here is an untested example based on the link you provided and radouxju's NDVI approach:
import numpy as np
dataset = r'C:\HY-DATA\HENTENKA\CSC\Data\Helsinki_masked_p188r018_7t20020529_z34__LV-FIN.tif'
# Outfile path
outpath = r'C:\path\to\...
The rasterio documentation for the mask function is not very clear here. The shape parameter must be an iterable of geometries, not a simple geometrie.
The documentation of the shape parameter (i.e. "The values must be a GeoJSON-like dict or an object that implements the Python geo interface protocol...") describes the nature of the values ...
You can iterate all over array, if that helps?
def reclassify(arr, rows, cols):
for r in range(rows):
for c in range(cols):
temp_value = arr[r][c]
if temp_value < 20:
temp_value = 4
elif temp_value >=80:
temp_value = 1
elif (temp_value >=50) and (...