05_load_raster¶
Following https://leafmap.org/notebooks/05_load_raster/
Uses the ipyleaflet plotting backend. The folium plotting backend does not support the add_raster function.
conda create -n leafmap -c conda-forge leafmap xarray_leaflet localtileserver netcdf4
conda activate leafmap
In [1]:
Copied!
import os
import leafmap.leafmap as leafmap
import os
import leafmap.leafmap as leafmap
load tif files from the example¶
Download sample raster datasets
More datasets can be downloaded from https://viewer.nationalmap.gov/basic/
In [2]:
Copied!
landsat = "landsat.tif"
dem = "dem.tif"
landsat = "landsat.tif"
dem = "dem.tif"
In [3]:
Copied!
url1 = "https://open.gishub.org/data/raster/landsat7.tif"
url2 = "https://open.gishub.org/data/raster/srtm90.tif"
satellite = leafmap.download_file(url1, landsat)
dem = leafmap.download_file(url2, dem)
url1 = "https://open.gishub.org/data/raster/landsat7.tif"
url2 = "https://open.gishub.org/data/raster/srtm90.tif"
satellite = leafmap.download_file(url1, landsat)
dem = leafmap.download_file(url2, dem)
landsat.tif already exists. Skip downloading. Set overwrite=True to overwrite. dem.tif already exists. Skip downloading. Set overwrite=True to overwrite.
In [4]:
Copied!
! gdalinfo landsat.tif
! gdalinfo landsat.tif
Driver: GTiff/GeoTIFF
Files: landsat.tif
Size is 2181, 1917
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Web mapping and visualisation."],
AREA["World between 85.06°S and 85.06°N."],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-13651650.000000000000000,4576290.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
TIFFTAG_RESOLUTIONUNIT=1 (unitless)
OVR_RESAMPLING_ALG=NEAREST
AREA_OR_POINT=Area
Image Structure Metadata:
LAYOUT=COG
COMPRESSION=DEFLATE
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-13651650.000, 4576290.000) (122d38' 5.49"W, 37d58'40.08"N)
Lower Left (-13651650.000, 4518780.000) (122d38' 5.49"W, 37d34'10.00"N)
Upper Right (-13586220.000, 4576290.000) (122d 2'49.53"W, 37d58'40.08"N)
Lower Right (-13586220.000, 4518780.000) (122d 2'49.53"W, 37d34'10.00"N)
Center (-13618935.000, 4547535.000) (122d20'27.51"W, 37d46'26.05"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
NoData Value=0
Overviews: 1091x959, 546x480
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
NoData Value=0
Overviews: 1091x959, 546x480
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
NoData Value=0
Overviews: 1091x959, 546x480
In [5]:
Copied!
m = leafmap.Map()
m.add_raster(dem, colormap="terrain", layer_name="DEM")
m.add_raster(landsat, bands=[1, 2, 3], layer_name="Landsat")
# m.add_raster(landsat, layer_name="Landsat")
m
m = leafmap.Map()
m.add_raster(dem, colormap="terrain", layer_name="DEM")
m.add_raster(landsat, bands=[1, 2, 3], layer_name="Landsat")
# m.add_raster(landsat, layer_name="Landsat")
m
Out[5]:
Map(center=[37.7736215, -122.34097449999999], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom…
load a netcdf file¶
To be reviewed..
In [8]:
Copied!
! ncdump -h ../qva/qva_det_Ruapehu_202501081200.nc
! ncdump -h ../qva/qva_det_Ruapehu_202501081200.nc
netcdf qva_det_Ruapehu_202501081200 {
dimensions:
time = UNLIMITED ; // (24 currently)
latitude = 361 ;
longitude = 441 ;
levels = 12 ;
variables:
float latitude(latitude) ;
latitude:_FillValue = NaNf ;
latitude:long_name = "latitude degrees north from the equator" ;
latitude:units = "degrees_north" ;
latitude:point_spacing = "even" ;
float longitude(longitude) ;
longitude:_FillValue = NaNf ;
longitude:long_name = "longitude degrees east from the greenwich meridian" ;
longitude:units = "degrees_east" ;
longitude:point_spacing = "even" ;
double levels(levels) ;
levels:_FillValue = NaN ;
levels:long_name = "Top of flight level layer" ;
levels:units = "feet" ;
levels:axis = "Z" ;
levels:positive = "up" ;
int64 time(time) ;
time:axis = "T" ;
time:long_name = "time" ;
time:units = "hours since 2025-01-08T13:00:00" ;
time:calendar = "proleptic_gregorian" ;
float concentration(time, levels, latitude, longitude) ;
concentration:_FillValue = NaNf ;
concentration:long_name = "Quantitative volcanic ash concentration from deterministic forecast" ;
concentration:units = "mg/m3" ;
// global attributes:
:title = "HYSPLIT Model Concentration Output" ;
:Conventions = "CF-1.5" ;
:volcano_name = "Ruapehu" ;
:eruption_lon_degrees = 175.57f ;
:eruption_lat_degrees = -39.28f ;
:eruption_vent_meters_msl = 2797.f ;
:eruption_height_meters_msl = 12500.f ;
:eruption_mass_eruption_rate_kgs = 87500. ;
:eruption_start_time = "2025-01-08T12:00:00" ;
:eruption_duration_s = 3600. ;
}
In [36]:
Copied!
ds = leafmap.read_netcdf('../qva/qva_det_Ruapehu_202501081200.nc')
ds.rio.write_crs("EPSG:4326", inplace=True)
da = ds['concentration'].isel(time=2, levels=7)
da = da.where(da > 1e-6)
da.plot()
ds = leafmap.read_netcdf('../qva/qva_det_Ruapehu_202501081200.nc')
ds.rio.write_crs("EPSG:4326", inplace=True)
da = ds['concentration'].isel(time=2, levels=7)
da = da.where(da > 1e-6)
da.plot()
Out[36]:
<matplotlib.collections.QuadMesh at 0x706ccb8bdd10>
In [37]:
Copied!
m = leafmap.Map()
m.add_raster(da, layer_name="QVA")
m
m = leafmap.Map()
m.add_raster(da, layer_name="QVA")
m
Out[37]:
Map(center=[-39.279999000000004, 175.570007], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom…
In [ ]:
Copied!