PRISMA and EMIT readers¶
Install package dependecies¶
EMIT requires the netcdf4 package to read the products. We will also use the pysolar package to convert the EMIT L1B radiances to reflectances. PRISMA requires the h5py package and we need scipy for the georeader.reflectance module.
pip install georeader-spaceml netcdf4 pysolar h5py scipy
Read images¶
import os
prisma_path = "prisma_database/PRS_L1_STD_OFFL_20230929102749_20230929102753_0001.he5"
emit_path = "emit_database/raw/EMIT_L1B_RAD_001_20230929T122534_2327208_039.nc"
os.path.exists(prisma_path),os.path.exists(emit_path)
from georeader.readers import emit
from georeader.readers import prisma
import matplotlib.pyplot as plt
from georeader import plot
from georeader import read
import matplotlib.pyplot as plt
from shapely.geometry import Point
from georeader.window_utils import polygon_to_crs
ei = emit.EMITImage(emit_path).to_crs()
ei
pi = prisma.PRISMA(prisma_path)
pi
# Load SWIR and NIR PRISMA images
pi.load_raw(swir_flag=True)
_ = pi.load_raw(swir_flag=False)
Read same pixel value¶
pixel_prisma = (45, 65)
# The raw values of prisma are col major
point_prisma = pi.lons[pixel_prisma[-1::-1]],pi.lats[pixel_prisma[-1::-1]]
ei_on_point = read.read_from_center_coords(ei, point_prisma,
shape=(50,50),
crs_center_coords="EPSG:4326").load()
ei_on_point
# Center pixel
pixel_value_emit = ei_on_point.values[:, 25, 25]
pixel_value_prisma = pi.ltoa_swir[pixel_prisma + (slice(None),)]
pixel_value_prisma_vnir = pi.ltoa_vnir[pixel_prisma + (slice(None),)]
plt.plot(ei.wavelengths, pixel_value_emit, label="EMIT")
plt.plot(pi.wavelength_swir[pixel_prisma[1],:], pixel_value_prisma/10, label="PRISMA SWIR")
plt.plot(pi.wavelength_vnir[pixel_prisma[1],:], pixel_value_prisma_vnir/10, label="PRISMA VNIR")
plt.ylabel(ei.units)
plt.xlabel("$\lambda$ (nm)")
plt.grid()
plt.legend()
Load RGB reflectance¶
rgb_prisma = pi.load_rgb(raw=False)
rgb_prisma
plot.show(rgb_prisma)
rgb_emit = ei.load_rgb(as_reflectance=True)
rgb_emit
fig, ax = plt.subplots(1,1)
plot.show(rgb_emit, ax=ax)
plot.add_shape_to_plot(ei.footprint(crs="EPSG:4326"), crs_plot=ei.crs, crs_shape="EPSG:4326", polygon_no_fill=True, ax=ax)
Plot RGB over the exact same area¶
import geopandas as gpd
metadata = gpd.GeoDataFrame({"geometry": [ei.footprint(crs="EPSG:4326"), pi.footprint(crs="EPSG:4326")],
# "geometry": [rgb_emit.valid_footprint(crs="EPSG:4326"), rgb_prisma.valid_footprint(crs="EPSG:4326")],
"satellite": ["EMIT", "PRISMA"],
"tile_date": [ei.time_coverage_start.isoformat(), pi.time_coverage_start.isoformat()]},
geometry="geometry",
crs="EPSG:4326")
metadata
metadata.explore()
pol_intersection = rgb_emit.valid_footprint(crs="EPSG:4326").intersection(rgb_prisma.valid_footprint(crs="EPSG:4326"))
rgb_emit_intersect = read.read_from_polygon(rgb_emit, pol_intersection, "EPSG:4326")
rgb_prisma_intersect = read.read_from_polygon(rgb_prisma, pol_intersection, "EPSG:4326")
# Show point where we extracted the profile in the plot
point_prisma_epsg = Point(*point_prisma)
point_emit_crs = polygon_to_crs(point_prisma_epsg, crs_polygon="EPSG:4326", dst_crs=rgb_emit_intersect.crs)
point_prisma_crs = polygon_to_crs(point_prisma_epsg, crs_polygon="EPSG:4326", dst_crs=rgb_prisma_intersect.crs)
fig, ax = plt.subplots(1,2,figsize=(16,10))
plot.show(rgb_emit_intersect, ax=ax[0], title="EMIT")
ax[0].scatter([point_emit_crs.coords[0][0]], [point_emit_crs.coords[0][1]], marker="x", color="red")
plot.add_shape_to_plot(pol_intersection, crs_plot=rgb_emit_intersect.crs, crs_shape="EPSG:4326", polygon_no_fill=True, ax=ax[0])
plot.show(rgb_prisma_intersect, ax=ax[1], title="PRISMA")
ax[1].scatter([point_prisma_crs.coords[0][0]], [point_prisma_crs.coords[0][1]], marker="x", color="red")
plot.add_shape_to_plot(pol_intersection, crs_plot=rgb_prisma_intersect.crs, crs_shape="EPSG:4326", polygon_no_fill=True, ax=ax[1])
Zoom in in the point¶
rgb_emit_point = read.read_from_center_coords(rgb_emit, point_prisma, shape=(50,50),
crs_center_coords="EPSG:4326")
rgb_prisma_point = read.read_from_center_coords(rgb_prisma, point_prisma,
shape=(100,100),
crs_center_coords="EPSG:4326")
# # Show point where we extracted the profile in the plot
# point_prisma_epsg = Point(*point_prisma)
# point_emit_crs = polygon_to_crs(point_prisma_epsg, crs_polygon="EPSG:4326", dst_crs=rgb_emit_intersect.crs)
# point_prisma_crs = polygon_to_crs(point_prisma_epsg, crs_polygon="EPSG:4326", dst_crs=rgb_prisma_intersect.crs)
fig, ax = plt.subplots(1,2,figsize=(16,10))
plot.show(rgb_emit_point, ax=ax[0], title="EMIT")
ax[0].scatter([point_emit_crs.coords[0][0]], [point_emit_crs.coords[0][1]], marker="x", color="red")
plot.show(rgb_prisma_point, ax=ax[1], title="PRISMA")
ax[1].scatter([point_prisma_crs.coords[0][0]], [point_prisma_crs.coords[0][1]], marker="x", color="red")
Licence¶
The georeader package is published under a GNU Lesser GPL v3 licence
If you find this work useful please cite:
@article{ruzicka_starcop_2023,
title = {Semantic segmentation of methane plumes with hyperspectral machine learning models},
volume = {13},
issn = {2045-2322},
url = {https://www.nature.com/articles/s41598-023-44918-6},
doi = {10.1038/s41598-023-44918-6},
number = {1},
journal = {Scientific Reports},
author = {Růžička, Vít and Mateo-Garcia, Gonzalo and Gómez-Chova, Luis and Vaughan, Anna, and Guanter, Luis and Markham, Andrew},
month = nov,
year = {2023},
pages = {19999},
}