89 image array viz
Visualizing in-memory raster datasets and image arrays
In [1]:
Copied!
# %pip install "leafmap[raster]"
# %pip install "leafmap[raster]"
In [2]:
Copied!
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]:
Copied!
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, 182MB/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]:
Copied!
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]:
Copied!
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]:
Copied!
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]:
Copied!
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]:
Copied!
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]:
Copied!
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]:
Copied!
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]:
Copied!
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]:
Copied!
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: elevationClassify the DEM into 2 elevation classes.
In [13]:
Copied!
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]:
Copied!
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]:
Copied!
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]: