Last active
July 3, 2025 19:07
-
-
Save previtus/d69019039a62d182df6827b2b037ccba to your computer and use it in GitHub Desktop.
EMIT snippets
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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"]) |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
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