Skip to main content Home Creators xjtulyc awesome-rosetta-skills rioxarray-remote-sensing
rioxarray-remote-sensing Read, reproject, resample, and analyze multispectral GeoTIFF imagery with rioxarray and rasterio; compute NDVI/EVI, mask clouds, mosaic scenes, and export Cloud-Optimized GeoTIFFs.
Jump to install Skills Marketplace Discover and explore AI skills built by the community.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
Show prompt details A direct command skips the review prompt. Inspect the source before running it.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill rioxarray-remote-sensingThe command stays on one line. Scroll horizontally to inspect it before copying.
Prefer a local copy? Download the files currently available to SkillsMP.
Download Zip Downloading... More from this repository
Related occupations SOC
Based on SOC occupation classification
name rioxarray-remote-sensing description Read, reproject, resample, and analyze multispectral GeoTIFF imagery with rioxarray and rasterio; compute NDVI/EVI, mask clouds, mosaic scenes, and export Cloud-Optimized GeoTIFFs.
tags ["remote-sensing","rioxarray","rasterio","ndvi","geotiff","cloud-optimized"] version 1.0.0 authors [{"name":"awesome-rosetta-skills contributors","github":"@xjtulyc"}] license MIT platforms ["claude-code","codex","gemini-cli","cursor"] dependencies {"python":["rioxarray>=0.15.0","rasterio>=1.3.9","earthpy>=0.9.4","numpy>=1.24.0","matplotlib>=3.7.0","geopandas>=0.14.0","shapely>=2.0.0"]} last_updated 2026-03-17 status stable
rioxarray Remote Sensing Skill
One-line summary : This Skill helps researchers read, reproject, resample,
and analyse multispectral satellite imagery using rioxarray/rasterio, suited
for NDVI/EVI vegetation analysis, cloud masking, mosaicking, and COG export.
When to Use This Skill
Use this Skill in the following scenarios:
When you need to read and write GeoTIFF files with full CRS/metadata preservation.
When your raster data needs reprojection or resampling to match a target grid.
When you need to compute vegetation indices (NDVI, EVI) from multispectral bands.
When you are using Landsat or Sentinel imagery with QA bit-mask cloud masking.
When you need to mosaic multiple scenes into a single seamless raster.
When you need to produce Cloud-Optimized GeoTIFFs (COGs) for web delivery.
When you need to clip rasters to polygon boundaries (watersheds, administrative units).
Trigger keywords : rioxarray, rasterio, NDVI, EVI, GeoTIFF, COG, CRS reprojection,
resampling, cloud masking, QA band, mosaicking, band math, clip to polygon.
Background & Key Concepts
Coordinate Reference Systems and Reprojection
A Coordinate Reference System (CRS) defines how 2-D map coordinates correspond to
locations on the Earth's surface. Satellite imagery is delivered in a variety of
projections — Landsat Collection 2 uses UTM zones, while global composites often use
WGS-84 geographic coordinates (EPSG:4326). Reprojection transforms pixel data from one
CRS to another by resampling values onto a new grid.
rioxarray wraps rasterio's GDAL-backed resampling. Available algorithms include:
Algorithm Use case Nearest neighbour Categorical data (land cover, QA bands) Bilinear Continuous data with moderate accuracy needs Cubic High-accuracy continuous data resampling Lanczos Downsampling with minimal aliasing
NDVI and EVI Formulae
The Normalized Difference Vegetation Index (NDVI) leverages the spectral contrast between
the Near-Infrared (NIR) and Red bands:
$$
\text{NDVI} = \frac{\rho_{NIR} - \rho_{Red}}{\rho_{NIR} + \rho_{Red}}
$$
Values range from -1 to +1; healthy vegetation typically falls between 0.3 and 0.8.
Sparse vegetation, bare soil, and water produce values closer to 0 or negative.
The Enhanced Vegetation Index (EVI) reduces atmospheric and soil-background noise:
$$
\text{EVI} = G \cdot \frac{\rho_{NIR} - \rho_{Red}}{\rho_{NIR} + C_1 \cdot \rho_{Red} - C_2 \cdot \rho_{Blue} + L}
$$
Where $G = 2.5$, $C_1 = 6$, $C_2 = 7.5$, $L = 1$ (standard MODIS coefficients).
Cloud Masking with QA Bits Landsat Collection 2 Quality Assessment (QA_PIXEL) uses bit-packed integers.
Bit 3 = cloud shadow, Bit 4 = cloud. Reading a specific bit:
$$
\text{bit}_k = \left\lfloor \frac{QA_PIXEL}{2^k} \right\rfloor \mod 2
$$
A pixel is considered cloud-contaminated when bit 3 OR bit 4 is set.
Comparison with Related Methods Method Best for Key assumption Limitation rioxarray + rasterio Python-native raster I/O GDAL-supported formats Python memory limits for very large mosaics GDAL CLI (gdalwarp) Batch reprojection, scripting Shell environment No Python integration Google Earth Engine Petabyte-scale cloud processing Internet + GEE account Proprietary platform xarray + rasterio NetCDF-heavy workflows Rectangular grids Less direct raster support
Environment Setup
Install Dependencies
python -m venv .venv
source .venv/bin/activate
pip install rioxarray>=0.15.0 rasterio>=1.3.9 earthpy>=0.9.4 \
numpy>=1.24.0 matplotlib>=3.7.0 geopandas>=0.14.0 shapely>=2.0.0
Verify Installation import rioxarray
import rasterio
import earthpy
import numpy as np
import matplotlib
import geopandas
print (f"rioxarray : {rioxarray.__version__} " )
print (f"rasterio : {rasterio.__version__} " )
print (f"numpy : {np.__version__} " )
print (f"geopandas : {geopandas.__version__} " )
Core Workflow
Step 1: Reading a GeoTIFF and CRS Inspection rioxarray's open_rasterio returns an xarray.DataArray with a .rio accessor that
exposes all CRS and spatial metadata.
import rioxarray as rxr
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import rasterio
from rasterio.transform import from_bounds
from rasterio.crs import CRS
np.random.seed(42 )
n = 64
synthetic = np.random.randint(500 , 8000 , (4 , n, n), dtype=np.int16)
synthetic[3 , 10 :40 , 10 :40 ] = np.random.randint(5000 , 9000 , (30 , 30 ))
synthetic[2 , 10 :40 , 10 :40 ] = np.random.randint(500 , 2000 , (30 , 30 ))
transform = from_bounds(500000 , 5700000 , 501920 , 5701920 , n, n)
crs = CRS.from_epsg(32632 )
scene_path = Path("synthetic_landsat_b2345.tif" )
with rasterio.open (
scene_path, "w" , driver="GTiff" ,
height=n, width=n, count=4 , dtype="int16" ,
crs=crs, transform=transform,
) as dst:
dst.write(synthetic)
dst.update_tags(1 , name="Blue" )
dst.update_tags(2 , name="Green" )
dst.update_tags(3 , name="Red" )
dst.update_tags(4 , name="NIR" )
print (f"Wrote synthetic scene: {scene_path} " )
da = rxr.open_rasterio(scene_path, masked=True )
print (f"Shape : {da.shape} " )
print (f"CRS : {da.rio.crs} " )
print (f"Resolution : {da.rio.resolution()} " )
print (f"Bounds : {da.rio.bounds()} " )
print (f"NoData value : {da.rio.nodata} " )
band: integer band index starting at 1y / x: spatial coordinate dimensions in the native CRS unitsRecommended: floating-point or 16-bit integer reflectance values, scale factor noted
Step 2: Reprojection and Resampling import rioxarray as rxr
from rasterio.enums import Resampling
da = rxr.open_rasterio("synthetic_landsat_b2345.tif" , masked=True )
da_wgs84 = da.rio.reproject("EPSG:4326" , resampling=Resampling.bilinear)
print (f"Reprojected CRS : {da_wgs84.rio.crs} " )
print (f"Reprojected bounds : {da_wgs84.rio.bounds()} " )
da_resampled = da.rio.reproject(
"EPSG:4326" ,
resolution=0.0003 ,
resampling=Resampling.cubic,
)
print (f"Custom-resolution shape: {da_resampled.shape} " )
reference = rxr.open_rasterio("synthetic_landsat_b2345.tif" , masked=True )
da_matched = da_wgs84.rio.reproject_match(reference)
print (f"Matched shape: {da_matched.shape} " )
Parameter Meaning Recommended value Notes resamplingPixel value interpolation method bilinear for continuousUse nearest for categorical layers resolutionTarget pixel size in CRS units Match target dataset Degrees for geographic CRS, metres for projected nodataFill value for areas outside original extent np.nanSet via da.rio.write_nodata() before reprojection
Step 3: NDVI and EVI Calculation import rioxarray as rxr
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore" )
da = rxr.open_rasterio("synthetic_landsat_b2345.tif" , masked=True ).astype(float )
blue = da.sel(band=1 ) / 10000.0
green = da.sel(band=2 ) / 10000.0
red = da.sel(band=3 ) / 10000.0
nir = da.sel(band=4 ) / 10000.0
ndvi = (nir - red) / (nir + red)
ndvi = ndvi.where(np.isfinite(ndvi))
G, C1, C2, L = 2.5 , 6.0 , 7.5 , 1.0
evi_denom = nir + C1 * red - C2 * blue + L
evi = G * (nir - red) / evi_denom
evi = evi.clip(-1.0 , 1.0 )
print (f"NDVI min/mean/max : {float (ndvi.min ()):.3 f} / "
f"{float (ndvi.mean()):.3 f} / {float (ndvi.max ()):.3 f} " )
print (f"EVI min/mean/max : {float (evi.min ()):.3 f} / "
f"{float (evi.mean()):.3 f} / {float (evi.max ()):.3 f} " )
fig, axes = plt.subplots(1 , 2 , figsize=(12 , 5 ))
im0 = axes[0 ].imshow(ndvi.values, cmap="RdYlGn" , vmin=-0.2 , vmax=0.9 )
axes[0 ].set_title("NDVI" )
axes[0 ].axis("off" )
plt.colorbar(im0, ax=axes[0 ], fraction=0.046 , pad=0.04 )
im1 = axes[1 ].imshow(evi.values, cmap="RdYlGn" , vmin=-0.2 , vmax=0.9 )
axes[1 ].set_title("EVI" )
axes[1 ].axis("off" )
plt.colorbar(im1, ax=axes[1 ], fraction=0.046 , pad=0.04 )
plt.suptitle("Vegetation Indices (synthetic Landsat-8 scene)" )
plt.tight_layout()
plt.savefig("ndvi_evi_comparison.png" , dpi=150 , bbox_inches="tight" )
plt.show()
print ("Saved: ndvi_evi_comparison.png" )
print ("\nVegetation cover fractions:" )
for label, lo, hi in [("Water/bare" , -1.0 , 0.1 ),
("Sparse veg" , 0.1 , 0.3 ),
("Moderate" , 0.3 , 0.6 ),
("Dense veg" , 0.6 , 1.0 )]:
frac = float (((ndvi >= lo) & (ndvi < hi)).sum ()) / ndvi.size
print (f" {label:<15 } : {frac*100 :5.1 f} %" )
NDVI > 0.6 : dense, healthy vegetation canopyNDVI 0.2–0.6 : moderate/sparse vegetationNDVI < 0.1 : bare soil, rock, urban areas, or water
Advanced Usage
Cloud Masking with QA Band (Bit Manipulation) Landsat Collection 2 QA_PIXEL band encodes quality flags as bit fields.
Bits 3 (cloud shadow) and 4 (cloud) indicate unreliable pixels.
import rioxarray as rxr
import numpy as np
import rasterio
from rasterio.transform import from_bounds
from rasterio.crs import CRS
from pathlib import Path
np.random.seed(0 )
n = 64
qa_data = np.zeros((1 , n, n), dtype=np.uint16)
cloud_mask_idx = np.random.choice(n * n, size=200 , replace=False )
rows, cols = np.unravel_index(cloud_mask_idx, (n, n))
qa_data[0 , rows, cols] |= (1 << 4 )
qa_data[0 , 5 :10 , 5 :10 ] |= (1 << 3 )
transform = from_bounds(500000 , 5700000 , 501920 , 5701920 , n, n)
qa_path = Path("synthetic_qa_pixel.tif" )
with rasterio.open (
qa_path, "w" , driver="GTiff" ,
height=n, width=n, count=1 , dtype="uint16" ,
crs=CRS.from_epsg(32632 ), transform=transform,
) as dst:
dst.write(qa_data)
qa = rxr.open_rasterio(qa_path, masked=False ).squeeze()
qa_np = qa.values.astype(np.uint16)
cloud_bit = 4
cloud_shadow_bit = 3
cloud_flag = (qa_np >> cloud_bit) & 1
cloud_shadow_flag = (qa_np >> cloud_shadow_bit) & 1
bad_pixel_mask = (cloud_flag | cloud_shadow_flag).astype(bool )
print (f"Flagged pixels: {bad_pixel_mask.sum ()} / {bad_pixel_mask.size} "
f"({bad_pixel_mask.mean()*100 :.1 f} %)" )
da = rxr.open_rasterio("synthetic_landsat_b2345.tif" , masked=True ).astype(float )
ndvi = (da.sel(band=4 ) - da.sel(band=3 )) / (da.sel(band=4 ) + da.sel(band=3 ))
ndvi_masked = ndvi.where(~bad_pixel_mask)
print (f"NDVI valid pixels after cloud masking: "
f"{int (np.isfinite(ndvi_masked.values).sum ())} " )
Clip Raster to Polygon Boundary import rioxarray as rxr
import geopandas as gpd
from shapely.geometry import box, mapping
import numpy as np
aoi = gpd.GeoDataFrame(
geometry=[box(500300 , 5700300 , 501600 , 5701600 )],
crs="EPSG:32632" ,
)
da = rxr.open_rasterio("synthetic_landsat_b2345.tif" , masked=True )
aoi = aoi.to_crs(da.rio.crs)
da_clipped = da.rio.clip(aoi.geometry.apply(mapping), aoi.crs, drop=True )
print (f"Original shape : {da.shape} " )
print (f"Clipped shape : {da_clipped.shape} " )
print (f"Clipped bounds : {da_clipped.rio.bounds()} " )
Mosaicking Multiple Scenes When a study area spans multiple raster tiles, merge them into a single seamless image.
import rioxarray as rxr
from rioxarray.merge import merge_arrays
import rasterio
import numpy as np
from rasterio.transform import from_bounds
from rasterio.crs import CRS
from pathlib import Path
def make_tile (path, x_min, x_max, n=64 ):
np.random.seed(int (x_min))
data = np.random.randint(500 , 8000 , (4 , n, n), dtype=np.int16)
transform = from_bounds(x_min, 5700000 , x_max, 5701920 , n, n)
with rasterio.open (
path, "w" , driver="GTiff" ,
height=n, width=n, count=4 , dtype="int16" ,
crs=CRS.from_epsg(32632 ), transform=transform,
) as dst:
dst.write(data)
make_tile("tile_west.tif" , 500000 , 501920 )
make_tile("tile_east.tif" , 501920 , 503840 )
tile_w = rxr.open_rasterio("tile_west.tif" , masked=True )
tile_e = rxr.open_rasterio("tile_east.tif" , masked=True )
mosaic = merge_arrays([tile_w, tile_e])
print (f"West tile shape : {tile_w.shape} " )
print (f"East tile shape : {tile_e.shape} " )
print (f"Mosaic shape : {mosaic.shape} " )
print (f"Mosaic bounds : {mosaic.rio.bounds()} " )
Export as Cloud-Optimized GeoTIFF (COG) COGs store overviews and use tiled storage for efficient HTTP range-request access.
import rioxarray as rxr
import numpy as np
da = rxr.open_rasterio("synthetic_landsat_b2345.tif" , masked=True )
nir = da.sel(band=4 ).astype("float32" )
red = da.sel(band=3 ).astype("float32" )
ndvi = (nir - red) / (nir + red + 1e-9 )
ndvi = ndvi.expand_dims("band" )
ndvi.rio.to_raster(
"ndvi_cog.tif" ,
driver="GTiff" ,
dtype="float32" ,
compress="deflate" ,
predictor=2 ,
tiled=True ,
blockxsize=512 ,
blockysize=512 ,
copy_src_overwrite=True ,
)
print ("Exported Cloud-Optimized GeoTIFF: ndvi_cog.tif" )
import rasterio
with rasterio.open ("ndvi_cog.tif" ) as src:
print (f" Driver : {src.driver} " )
print (f" CRS : {src.crs} " )
print (f" Shape : {src.width} x {src.height} " )
print (f" Profile : {src.profile} " )
Performance Optimization import rioxarray as rxr
import dask
da_lazy = rxr.open_rasterio(
"large_scene.tif" ,
masked=True ,
chunks={"band" : 1 , "x" : 1024 , "y" : 1024 },
lock=False ,
)
print (f"Dask DataArray: {da_lazy} " )
print (f"Chunk sizes : {da_lazy.chunks} " )
da_subset = da_lazy.isel(x=slice (0 , 2048 ), y=slice (0 , 2048 )).compute()
print (f"Subset shape : {da_subset.shape} " )
Troubleshooting
Error: CRSError: Invalid projection: EPSG:XXXX Cause : EPSG code not found in local PROJ database; proj.db may be outdated.
conda install -c conda-forge proj>=9.0
pip install --force-reinstall rasterio
Error: NoDataError after clipping Cause : The clip polygon does not overlap the raster extent, or CRS mismatch.
import rioxarray as rxr
import geopandas as gpd
da = rxr.open_rasterio("my_raster.tif" )
gdf = gpd.read_file("my_polygon.shp" )
gdf_reproj = gdf.to_crs(da.rio.crs)
print (f"Raster CRS : {da.rio.crs} " )
print (f"Polygon CRS : {gdf_reproj.crs} " )
print (f"Raster bounds : {da.rio.bounds()} " )
print (f"Polygon bounds : {gdf_reproj.total_bounds} " )
Issue: NDVI values outside [-1, 1] Cause : Unscaled integer reflectance values (e.g., Landsat SR 0–10000 range)
passed directly into the NDVI formula.
Fix : Divide by the scale factor (10000 for Landsat C2 L2 SR) before computing
band math.
Version Compatibility Package Tested versions Known issues rioxarray 0.15.x, 0.16.x reproject_match signature changed in 0.14; use keyword argsrasterio 1.3.x, 1.4.x Windows wheels require GDAL ≥ 3.7; use conda on Windows GDAL (system) 3.7–3.9 COG creation requires GDAL ≥ 3.1 numpy 1.24–2.0 Boolean indexing behaviour stable across versions
External Resources
Official Documentation
Key Papers
Rouse et al. (1974). Monitoring Vegetation Systems in the Great Plains with ERTS. NASA.
Huete et al. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83, 195–213.
Cloud Optimized GeoTIFF specification: https://cogeo.org/
Tutorials
Data Sources
Examples
Example 1: Full NDVI Processing Pipeline for a Landsat Scene Scenario : Compute cloud-masked NDVI for a Landsat-8 Collection 2 Level-2 scene,
clip to an agricultural region of interest, and export as COG.
Input data : Landsat-8 Collection 2 Level-2 SR bands (B2–B5) + QA_PIXEL band.
The example below uses synthetic data to be fully runnable without downloading.
import numpy as np
import rasterio
import rioxarray as rxr
from rasterio.transform import from_bounds
from rasterio.crs import CRS
from rioxarray.merge import merge_arrays
import geopandas as gpd
from shapely.geometry import box, mapping
import matplotlib.pyplot as plt
from pathlib import Path
import warnings
warnings.filterwarnings("ignore" )
np.random.seed(123 )
n = 128
transform = from_bounds(500000 , 5700000 , 503840 , 5703840 , n, n)
crs = CRS.from_epsg(32632 )
reflectance = np.random.randint(200 , 4000 , (4 , n, n), dtype=np.int16)
reflectance[3 , 30 :80 , 30 :80 ] = np.random.randint(6000 , 9500 , (50 , 50 ))
reflectance[2 , 30 :80 , 30 :80 ] = np.random.randint(300 , 1200 , (50 , 50 ))
qa = np.zeros((1 , n, n), dtype=np.uint16)
cloud_y, cloud_x = np.random.choice(n, 300 ), np.random.choice(n, 300 )
qa[0 , cloud_y, cloud_x] |= (1 << 4 )
for fname, data, dtype in [
("ls8_sr.tif" , reflectance, "int16" ),
("ls8_qa.tif" , qa, "uint16" ),
]:
with rasterio.open (fname, "w" , driver="GTiff" ,
height=n, width=n,
count=data.shape[0 ], dtype=dtype,
crs=crs, transform=transform) as dst:
dst.write(data)
da = rxr.open_rasterio("ls8_sr.tif" , masked=True ).astype("float32" )
qa_da = rxr.open_rasterio("ls8_qa.tif" , masked=False ).squeeze()
qa_np = qa_da.values.astype(np.uint16)
cloud_flag = (qa_np >> 4 ) & 1
cloud_shadow_flag = (qa_np >> 3 ) & 1
bad = (cloud_flag | cloud_shadow_flag).astype(bool )
print (f"Cloud/shadow pixels: {bad.sum ()} ({bad.mean()*100 :.1 f} %)" )
nir = da.sel(band=4 ) / 10000.0
red = da.sel(band=3 ) / 10000.0
ndvi = (nir - red) / (nir + red)
ndvi_clean = ndvi.where(~bad)
aoi = gpd.GeoDataFrame(geometry=[box(500500 , 5700500 , 503000 , 5703000 )], crs=crs)
ndvi_aoi = ndvi_clean.expand_dims("band" ).rio.clip(
aoi.geometry.apply(mapping), aoi.crs, drop=True
).squeeze("band" )
fig, axes = plt.subplots(1 , 2 , figsize=(12 , 5 ))
axes[0 ].imshow(ndvi_clean.values, cmap="RdYlGn" , vmin=-0.1 , vmax=0.9 )
axes[0 ].set_title("NDVI (cloud-masked, full scene)" )
axes[0 ].axis("off" )
axes[1 ].imshow(ndvi_aoi.values, cmap="RdYlGn" , vmin=-0.1 , vmax=0.9 )
axes[1 ].set_title("NDVI (clipped to AOI)" )
axes[1 ].axis("off" )
plt.tight_layout()
plt.savefig("ndvi_pipeline_output.png" , dpi=150 )
plt.show()
ndvi_aoi.expand_dims("band" ).rio.to_raster(
"ndvi_aoi_cog.tif" ,
driver="GTiff" , dtype="float32" ,
compress="deflate" , tiled=True , blockxsize=256 , blockysize=256 ,
)
print (f"NDVI stats (AOI): min={float (ndvi_aoi.min ()):.3 f} , "
f"mean={float (ndvi_aoi.mean()):.3 f} , max={float (ndvi_aoi.max ()):.3 f} " )
print ("Exported: ndvi_aoi_cog.tif" )
Interpreting these results : The vegetated patch (rows 30–80, cols 30–80) should
show NDVI around 0.6–0.8. Cloud-masked pixels appear as NaN (white in the figure).
The COG is ready for web map tile services.
Example 2: Multi-Scene Mosaic and EVI Map Scenario : Combine two adjacent Landsat tiles and produce a gap-free EVI mosaic.
import numpy as np
import rasterio
import rioxarray as rxr
from rioxarray.merge import merge_arrays
from rasterio.transform import from_bounds
from rasterio.crs import CRS
import matplotlib.pyplot as plt
def make_scene (path, x_min, x_max, seed=0 ):
np.random.seed(seed)
n = 64
data = np.random.randint(200 , 7000 , (4 , n, n), dtype=np.int16)
data[3 , 10 :40 , 10 :40 ] = np.random.randint(5500 , 9000 , (30 , 30 ))
data[2 , 10 :40 , 10 :40 ] = np.random.randint(400 , 1500 , (30 , 30 ))
data[0 , 10 :40 , 10 :40 ] = np.random.randint(200 , 700 , (30 , 30 ))
t = from_bounds(x_min, 5700000 , x_max, 5701920 , n, n)
with rasterio.open (path, "w" , driver="GTiff" ,
height=n, width=n, count=4 , dtype="int16" ,
crs=CRS.from_epsg(32632 ), transform=t) as dst:
dst.write(data)
make_scene("scene_A.tif" , 500000 , 501920 , seed=7 )
make_scene("scene_B.tif" , 501920 , 503840 , seed=13 )
tiles = [rxr.open_rasterio(p, masked=True ).astype("float32" )
for p in ["scene_A.tif" , "scene_B.tif" ]]
mosaic = merge_arrays(tiles)
print (f"Mosaic shape: {mosaic.shape} , bounds: {mosaic.rio.bounds()} " )
blue = mosaic.sel(band=1 ) / 10000.0
red = mosaic.sel(band=3 ) / 10000.0
nir = mosaic.sel(band=4 ) / 10000.0
G, C1, C2, L = 2.5 , 6.0 , 7.5 , 1.0
evi = G * (nir - red) / (nir + C1 * red - C2 * blue + L)
evi = evi.clip(-1.0 , 1.0 )
fig, ax = plt.subplots(figsize=(10 , 4 ))
im = ax.imshow(evi.values, cmap="RdYlGn" , vmin=-0.1 , vmax=0.9 , aspect="auto" )
ax.set_title("EVI — Mosaicked Scene (2 tiles)" )
ax.axis("off" )
plt.colorbar(im, ax=ax, fraction=0.02 , pad=0.04 , label="EVI" )
plt.tight_layout()
plt.savefig("evi_mosaic.png" , dpi=150 )
plt.show()
print (f"EVI: min={float (evi.min ()):.3 f} , mean={float (evi.mean()):.3 f} , "
f"max={float (evi.max ()):.3 f} " )
Interpreting these results : The mosaic seamlessly joins both tiles.
EVI values in the simulated vegetation patches (top-left quadrant of each tile)
should exceed 0.4. EVI is preferred over NDVI in dense canopy conditions because
it saturates less in high-biomass areas.
Last updated: 2026-03-17 | Maintainer: @xjtulyc
Issues: GitHub Issues
