Skip to content

89 image array viz

image image image

Visualizing in-memory raster datasets and image arrays

1
# %pip install "leafmap[raster]"

1
2
3
4
import leafmap
import rasterio
import rioxarray
import xarray as xr

Download two sample raster datasets.

1
2
3
4
url1 = "https://opengeos.org/data/raster/landsat.tif"
url2 = "https://opengeos.org/data/raster/srtm90.tif"
satellite = leafmap.download_file(url1, "landsat.tif", overwrite=True)
dem = leafmap.download_file(url2, "srtm90.tif")

The Landsat image contains 3 bands: nir, red, and green. Let's calculate NDVI using the nir and red bands.

1
2
3
4
dataset = rasterio.open(satellite)
nir = dataset.read(4).astype(float)
red = dataset.read(1).astype(float)
ndvi = (nir - red) / (nir + red)

Create an in-memory raster dataset from the NDVI array and use the projection and extent of the Landsat image.

1
ndvi_image = leafmap.array_to_image(ndvi, source=satellite)

Visualize the Landsat image and the NDVI image on the same map.

1
2
3
4
m = leafmap.Map()
m.add_raster(satellite, indexes=[4, 1, 2], vmin=0, vmax=120, layer_name="Landsat 7")
m.add_raster(ndvi_image, colormap="Greens", layer_name="NDVI")
m

You can also specify the image metadata (e.g., cellsize, crs, and transform) when creating the in-memory raster dataset.

First, check the metadata of the origina image.

1
dataset.profile

Check the crs of the original image.

1
dataset.crs

Check the transform of the original image.

1
dataset.transform

Create an in-memory raster dataset from the NDVI array and specify the cellsize, crs, and transform.

1
2
3
4
transform = (30.0, 0.0, -13651650.0, 0.0, -30.0, 4576290.0)
ndvi_image = leafmap.array_to_image(
    ndvi, cellsize=30, crs="EPSG:3857", transform=transform
)

Add the NDVI image to the map.

1
2
3
4
m = leafmap.Map()
m.add_raster(satellite, indexes=[4, 1, 2], vmin=0, vmax=120, layer_name="Landsat 7")
m.add_raster(ndvi_image, colormap="Greens", layer_name="NDVI")
m

Use rioxarray to read raster datasets into xarray DataArrays.

1
2
ds = rioxarray.open_rasterio(dem)
ds

Classify the DEM into 2 elevation classes.

1
2
array = ds.sel(band=1)
masked_array = xr.where(array < 2000, 0, 1)

Visualize the DEM and the elevation class image on the same map.

1
2
3
4
m = leafmap.Map()
m.add_raster(dem, colormap="terrain", layer_name="DEM")
m.add_raster(masked_array, colormap="coolwarm", layer_name="Classified DEM")
m

Add a split map.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
m = leafmap.Map(center=[37.6, -119], zoom=9)
m.split_map(
    dem,
    masked_array,
    left_args={
        "layer_name": "DEM",
        "colormap": "terrain",
    },
    right_args={
        "layer_name": "Classified DEM",
        "colormap": "coolwarm",
    },
)
m