89 image array viz

Visualizing in-memory raster datasets and image arrays

In [1]:

# %pip install "leafmap[raster]"

# %pip install "leafmap[raster]"

In [2]:

import leafmap
import rasterio
import rioxarray
import xarray as xr

import leafmap import rasterio import rioxarray import xarray as xr

Download two sample raster datasets.

In [3]:

url1 = "https://open.gishub.org/data/raster/landsat.tif"
url2 = "https://open.gishub.org/data/raster/srtm90.tif"
satellite = leafmap.download_file(url1, "landsat.tif", overwrite=True)
dem = leafmap.download_file(url2, "srtm90.tif")

url1 = "https://open.gishub.org/data/raster/landsat.tif" url2 = "https://open.gishub.org/data/raster/srtm90.tif" satellite = leafmap.download_file(url1, "landsat.tif", overwrite=True) dem = leafmap.download_file(url2, "srtm90.tif")

Downloading...
From: https://open.gishub.org/data/raster/landsat.tif
To: /home/runner/work/leafmap/leafmap/docs/notebooks/landsat.tif

  0%|          | 0.00/10.1M [00:00<?, ?B/s]

100%|██████████| 10.1M/10.1M [00:00<00:00, 293MB/s]

srtm90.tif already exists. Skip downloading. Set overwrite=True to overwrite.

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

In [4]:

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

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.

In [5]:

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

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

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

In [6]:

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

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

Out[6]:

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.

In [7]:

dataset.profile

dataset.profile

Out[7]:

{'driver': 'GTiff', 'dtype': 'uint8', 'nodata': 0.0, 'width': 2127, 'height': 1564, 'count': 4, 'crs': CRS.from_wkt('PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]'), 'transform': Affine(180.0, 0.0, -13442580.0,
       0.0, -180.0, 4670100.0), 'blockxsize': 512, 'blockysize': 512, 'tiled': True, 'compress': 'deflate', 'interleave': 'pixel'}

Check the crs of the original image.

In [8]:

dataset.crs

dataset.crs

Out[8]:

CRS.from_wkt('PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]')

Check the transform of the original image.

In [9]:

dataset.transform

dataset.transform

Out[9]:

Affine(180.0, 0.0, -13442580.0,
       0.0, -180.0, 4670100.0)

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

In [10]:

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
)

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.

In [11]:

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

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

Out[11]:

Use rioxarray to read raster datasets into xarray DataArrays.

In [12]:

ds = rioxarray.open_rasterio(dem)
ds

ds = rioxarray.open_rasterio(dem) ds

Out[12]:

<xarray.DataArray (band: 1, y: 2465, x: 4269)> Size: 21MB
[10523085 values with dtype=int16]
Coordinates:
  * band         (band) int64 8B 1
  * x            (x) float64 34kB -120.8 -120.8 -120.8 ... -117.3 -117.3 -117.3
  * y            (y) float64 20kB 38.63 38.63 38.62 38.62 ... 36.64 36.64 36.63
    spatial_ref  int64 8B 0
Attributes:
    OVR_RESAMPLING_ALG:  NEAREST
    AREA_OR_POINT:       Area
    scale_factor:        1.0
    add_offset:          0.0
    long_name:           elevation

xarray.DataArray

• band: 1
• y: 2465
• x: 4269

• ...

  [10523085 values with dtype=int16]

• Coordinates: (4)
  • band
    (band)
    int64
    1

    array([1])

  • x
    (x)
    float64
    -120.8 -120.8 ... -117.3 -117.3

    array([-120.755523, -120.754715, -120.753906, ..., -117.306531, -117.305723,
           -117.304914])

  • y
    (y)
    float64
    38.63 38.63 38.62 ... 36.64 36.63

    array([38.626519, 38.625711, 38.624902, ..., 36.636032, 36.635224, 36.634415])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -120.75592734926073 0.0008084837557075694 0.0 38.6269234341147 0.0 -0.0008084837557075694

    array(0)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index([1], dtype='int64', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([-120.75552310738287, -120.75471462362717, -120.75390613987146,
           -120.75309765611574, -120.75228917236004, -120.75148068860433,
           -120.75067220484863, -120.74986372109292, -120.74905523733722,
            -120.7482467535815,
           ...
           -117.31219079182434, -117.31138230806863, -117.31057382431291,
            -117.3097653405572,  -117.3089568568015,  -117.3081483730458,
           -117.30733988929009, -117.30653140553439, -117.30572292177867,
           -117.30491443802296],
          dtype='float64', name='x', length=4269))

  • y
    PandasIndex

    PandasIndex(Index([ 38.62651919223685, 38.625710708481144,  38.62490222472543,
            38.62409374096973,  38.62328525721402,  38.62247677345831,
           38.621668289702605,  38.62085980594689,  38.62005132219119,
            38.61924283843548,
           ...
            36.64169157197477, 36.640883088219056,  36.64007460446335,
           36.639266120707646, 36.638457636951934,  36.63764915319623,
           36.636840669440524,  36.63603218568481,  36.63522370192911,
           36.634415218173395],
          dtype='float64', name='y', length=2465))

• Attributes: (5)

  OVR_RESAMPLING_ALG :

  NEAREST

  AREA_OR_POINT :

  Area

  scale_factor :

  1.0

  add_offset :

  0.0

  long_name :

  elevation

Classify the DEM into 2 elevation classes.

In [13]:

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

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.

In [14]:

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

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

Out[14]:

Add a split map.

In [15]:

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

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

Out[15]: