Skip to content

Instantly share code, notes, and snippets.

@previtus
Last active July 3, 2025 19:07
Show Gist options
  • Save previtus/d69019039a62d182df6827b2b037ccba to your computer and use it in GitHub Desktop.
Save previtus/d69019039a62d182df6827b2b037ccba to your computer and use it in GitHub Desktop.
EMIT snippets
import matplotlib.pyplot as plt
import rasterio as rio
from georeader.geotensor import GeoTensor
from georeader.save import save_cog
import georeader
from georeader import read
from georeader.rasterio_reader import RasterioReader
from georeader.readers import emit
import numpy as np
import os
folder = "/Users/ruzicka/Downloads/DATA/OxHyper/_aux_data_for_ox_repos/"
label_l2b = folder+"EMIT_L2B_CH4PLM_001_20230612T162103_000320.tif"
data_l1b = folder+"EMIT_L1B_RAD_001_20230612T162103_2316311_006.nc"
def load_label_in_coord_of_emit(emit_path, label_path, in_utm_coords=True):
ei = emit.EMITImage(emit_path)
if in_utm_coords: ei = ei.to_crs(georeader.get_utm_epsg(ei.footprint("EPSG:4326")))
label_geo = RasterioReader(label_path)
label_reprojected = read.read_reproject_like(label_geo, ei, resampling = rio.warp.Resampling.nearest )
return label_reprojected.values, label_reprojected
def load_emit_data_in_location_of_label(emit_path, label_path, in_utm_coords=True):
ei = emit.EMITImage(emit_path)
if in_utm_coords: ei = ei.to_crs(georeader.get_utm_epsg(ei.footprint("EPSG:4326")))
label_geo = RasterioReader(label_path)
emit_data_reprojected = read.read_reproject_like(ei, label_geo, resampling = rio.warp.Resampling.nearest )
return emit_data_reprojected.values, emit_data_reprojected
def extract_tile_from_emit_data(emit_path, label_path, tile_size = 128, in_utm_coords=True):
ei = emit.EMITImage(emit_path)
if in_utm_coords: ei = ei.to_crs(georeader.get_utm_epsg(ei.footprint("EPSG:4326")))
label_geo = RasterioReader(label_path)
center_coords = label_geo.footprint("EPSG:4326").centroid.coords[0]
emit_data_tile = read.read_from_center_coords(ei, center_coords, shape=(tile_size, tile_size), crs_center_coords="EPSG:4326")
emit_data_tile = emit_data_tile.load(boundless = True, as_reflectance = False)
label_in_tile = read.read_from_center_coords(label_geo, center_coords, shape=(tile_size, tile_size), crs_center_coords="EPSG:4326")
return emit_data_tile.values, emit_data_tile, label_in_tile.values, label_in_tile
# First, let's show the label in the original EMIT tile's coords:
label_in_emit_coords, geo_obj = load_label_in_coord_of_emit(data_l1b, label_l2b)
print(label_in_emit_coords.shape)
plt.imshow(label_in_emit_coords[0])
plt.show()
save_cog(geo_obj, "label_in_emit_coords.tif", descriptions=["label"])
# Second, let's load EMIT data only in the location of the label:
emit_data_in_label, geo_obj = load_emit_data_in_location_of_label(data_l1b, label_l2b)
print(emit_data_in_label.shape) # (285, 36, 62)
rgb_indices = [35, 23, 11]
rgb_data = emit_data_in_label[rgb_indices, :,:]
plt.imshow(np.transpose(rgb_data / np.nanmax(rgb_data), (1, 2, 0)))
plt.show()
save_cog(geo_obj, "emit_data_in_label.tif") # < might want to add the descriptions
# Thirdly, we might want to load a full tile centered on this area
tile_data, tile_geo_obj, label_data, label_geo_obj = extract_tile_from_emit_data(data_l1b, label_l2b, tile_size = 128)
print(tile_data.shape) # (285, 36, 62)
rgb_data = tile_data[rgb_indices, :,:]
plt.imshow(np.transpose(rgb_data / np.nanmax(rgb_data), (1, 2, 0)))
plt.show()
plt.imshow(label_data[0])
plt.show()
save_cog(tile_geo_obj, "tile_data.tif") # < might want to add the descriptions
save_cog(label_geo_obj, "tile_label.tif", descriptions=["label"])
import os, glob
import rasterio
from tqdm import tqdm
if __name__ == "__main__":
plume_files_p = "/media/vitek/WD_BLACK/EMIT_raws/RAW_REAL_EVENTS/all_CH4PLM_Jan24"
plume_files = glob.glob(os.path.join(plume_files_p,"*.tif"))
print("we got", len(plume_files), "plume files")
# we got 862 plume files
complexes = []
sourcescenes = []
# now do we have dualities in between these?
for path in tqdm(plume_files):
with rasterio.open(path) as src:
# data = src.read()
t = src.tags()
complexes.append(t['Plume_Complex'])
sourcescenes.append(t['Source_Scenes'])
#print(t['Plume_Complex'])
print("Plume_Complex / in these we get", len(complexes), "-> unique:", len(list(set(complexes))))
# Plume_Complex / in these we get 862 -> unique: 862
# ^ all unique!
print("Source_Scenes / in these we get", len(sourcescenes), "-> unique:", len(list(set(sourcescenes))))
# Source_Scenes / in these we get 862 -> unique: 521
"""
Follows beer lambert law
- Garcia https://amt.copernicus.org/articles/15/1657/2022/
Xhat = X * e ^ (- AMF * \sigma_{CH_4} * mf)
(in our case) t * mf
- Foote https://ieeexplore.ieee.org/document/9034492
L(alpha, s) = L_bg * exp( k * alpha * s )
"""
def sim(X,t,y):
Xhat = np.zeros_like(X)
k = 1.
for band_i in range(len(X)):
Xhat[ band_i ] = X[ band_i ] * (np.exp( k * y * t[ band_i ] ) )
return Xhat
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment