Merge pull request #7 from emit-sds/develop

Merge develop into main for v1.3.0

CHANGELOG.md

@@ -4,10 +4,17 @@ All notable changes to this project will be documented in this file. Dates are d
 
 Generated by [`auto-changelog`](https://github.com/CookPete/auto-changelog).
 
+#### [v1.3.0](https://github.com/emit-sds/emit-sds-l3/compare/v1.2.0...v1.3.0)
+
+> 25 July 2024
+
+- add some v002 delivery scripts [`#6`](https://github.com/emit-sds/emit-sds-l3/pull/6)
+
 #### [v1.2.0](https://github.com/emit-sds/emit-sds-l3/compare/v1.1.0...v1.2.0)
 
 > 3 July 2024
 
+- Merge develop into main for v1.2.0 [`#5`](https://github.com/emit-sds/emit-sds-l3/pull/5)
 - add citations [`#4`](https://github.com/emit-sds/emit-sds-l3/pull/4)
 - Serial glt [`#3`](https://github.com/emit-sds/emit-sds-l3/pull/3)
 - add multithreading capability [`#2`](https://github.com/emit-sds/emit-sds-l3/pull/2)

asa_quicklook.py

@@ -0,0 +1,63 @@
+import argparse
+from mpl_toolkits.basemap import Basemap
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import matplotlib.gridspec as gridspec
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import matplotlib
+import glob
+from spectral.io import envi
+from osgeo import gdal
+
+
+deliver_band_names = [
+    "Calcite",
+    "Chlorite",
+    "Dolomite",
+    "Goethite",
+    "Gypsum",
+    "Hematite",
+    "Illite+Muscovite",
+    "Kaolinite",
+    "Montmorillonite",
+    "Vermiculite",
+]
+
+def main():
+
+    parser = argparse.ArgumentParser(description="Translate to Rrs. and/or apply masks")
+    parser.add_argument('input_file', type=str, metavar='as a nc')
+    parser.add_argument('output_file', type=str, metavar='output file to write')
+    parser.add_argument('--dpi', type=str, default='200')
+    args = parser.parse_args()
+
+
+    matplotlib.rc('font', **{'size'   : 22})
+    fig = plt.figure(figsize=(25, 20), edgecolor='w')
+    gs = gridspec.GridSpec(5, 2)
+
+    dat = envi.open(args.input_file + '.hdr').open_memmap(interleave='bip').copy()
+    trans = gdal.Open(args.input_file).GetGeoTransform()
+    band_names = envi.open(args.input_file + '.hdr').metadata['band names']
+
+    for n in range(len(deliver_band_names)):
+        ax = plt.subplot(gs[n % 5,int(np.floor(n/5))])
+
+        map = Basemap(llcrnrlon=trans[0], llcrnrlat=trans[3]+dat.shape[0]*trans[5], urcrnrlon=trans[0]+dat.shape[1]*trans[1], urcrnrlat=trans[3],
+                    lat_0=-0.01, lon_0=15.0, projection="eqdc")
+        map.drawcoastlines()
+
+
+        dat[dat == -9999] = np.nan
+        ldat = np.sum(dat[...,np.array([deliver_band_names[n] in x for x in band_names])],axis=-1)
+        maxval =max(np.nanpercentile(ldat,95), 0.05)
+        im = map.imshow(ldat[::-1,...],vmin=0,vmax=maxval)
+        plt.title(deliver_band_names[n] + ' Spectral Abundance')
+        cax1 = make_axes_locatable(ax).append_axes("right", size="2%", pad=0.02)
+        fig.colorbar(im, cax=cax1);
+
+    plt.savefig(args.output_file,dpi=200,bbox_inches='tight')
+
+if __name__ == "__main__":
+    main()

merge_for_delivery.py

@@ -0,0 +1,114 @@
+
+
+
+from netCDF4 import Dataset
+import argparse
+import numpy as np
+import os
+from datetime import datetime
+import pandas as pd
+from osgeo import osr
+import subprocess
+
+
+NODATA = -9999
+
+
+def main():
+    parser = argparse.ArgumentParser(description='netcdf conversion')
+    parser.add_argument('baseline_file', type=str)
+    parser.add_argument('unc_lower_file', type=str)
+    parser.add_argument('unc_upper_file', type=str)
+    parser.add_argument('output_file', type=str)
+    args = parser.parse_args()
+
+
+    baseline = Dataset(args.baseline_file, 'r')
+    uq_l = Dataset(args.unc_lower_file, 'r')
+    uq_u = Dataset(args.unc_upper_file, 'r')
+
+    output_dataset = Dataset(args.output_file, 'w', close=False, format='NETCDF4')
+    output_dataset.setncatts(baseline.__dict__)
+    for dname, dim in baseline.dimensions.items():
+        output_dataset.createDimension(dname, len(dim) if not dim.isunlimited() else None)
+
+
+    output_dataset.title = "EMIT L3 Aggregated Mineral Spectral Abundance 0.5 Deg. V002"
+    output_dataset.date_created = "2024-07-24T00:00:00Z"
+    output_dataset.license = "Freely Distributed"
+    output_dataset.creator_url = "https://www.jpl.nasa.gov"
+    output_dataset.summary = "The Earth Surface Mineral Dust Source Investigation (EMIT) is an Earth Ventures-Instrument (EVI-4) \
+Mission that maps the surface mineralogy of arid dust source regions via imaging spectroscopy in the visible and \
+short-wave infrared (VSWIR). Installed on the International Space Station (ISS), the EMIT instrument is a Dyson \
+imaging spectrometer that uses contiguous spectroscopic measurements from 410 to 2450 nm to resolve absoprtion \
+features of iron oxides, clays, sulfates, carbonates, and other dust-forming minerals. During its one-year mission, \
+EMIT will observe the sunlit Earth's dust source regions that occur within +/-52° latitude and produce maps of the \
+source regions that can be used to improve forecasts of the role of mineral dust in the radiative forcing \
+(warming or cooling) of the atmosphere.\n"
+    output_dataset.summary = output_dataset.summary + \
+    f"  This collection contains L3 Aggregated Mineral Spectral Abundance (ASA), at 0.5 degree resolution, \
+for use in Earth System Models.  ASA has been masked in areas with high vegetation, water, cloud, or urban cover. \
+The primary driver of uncertainty in the mass fractions is the grainsize used in the model.  As such, the uncertainty \
+is represented through the 2.5 and 97.5 percentiles of the grainsize distribution, designated as UQ_low_grainsize and UQ_high_grainsize."
+
+
+    for name, variable in baseline.variables.items():
+      if name not in ['latitude','longitude']:
+        new_var = output_dataset.createVariable(name, variable.datatype, variable.dimensions)
+        new_var.setncatts(variable.__dict__)
+        new_var.long_name = name + "_Baseline"
+        new_var[:] = variable[:]
+
+    for name, variable in uq_u.variables.items():
+      if name not in ['latitude','longitude']:
+        new_var = output_dataset.createVariable(name + "_UQ_high_grainsize", variable.datatype, variable.dimensions)
+        new_var.setncatts(variable.__dict__)
+        new_var.long_name = name + "_UQ_high_grainsize"
+        new_var[:] = variable[:]
+
+    for name, variable in uq_l.variables.items():
+      if name not in ['latitude','longitude']:
+        new_var = output_dataset.createVariable(name + "_UQ_low_grainsize", variable.datatype, variable.dimensions)
+        new_var.setncatts(variable.__dict__)
+        new_var.long_name = name + "_UQ_low_grainsize"
+        new_var[:] = variable[:]
+
+    new_var = output_dataset.createVariable('lat', baseline.variables['latitude'].datatype, baseline.variables['latitude'].dimensions)
+    new_var.standard_name = "latitude"
+    new_var.setncatts(baseline.variables['latitude'].__dict__)
+    new_var[:] = baseline.variables['latitude'][:]
+
+    new_var = output_dataset.createVariable('lon', baseline.variables['longitude'].datatype, baseline.variables['longitude'].dimensions)
+    new_var.standard_name = "longitude"
+    new_var.setncatts(baseline.variables['longitude'].__dict__)
+    new_var[:] = baseline.variables['longitude'][:]
+
+    # Add grid mapping variable if doesn't exist
+    if 'latitude_longitude' not in output_dataset.dimensions:
+        grid_mapping = output_dataset.createVariable('latitude_longitude', 'i4')
+
+
+        lat = np.sort(output_dataset.variables['lat'])
+        lon = np.sort(output_dataset.variables['lon'])
+        dlat = lat[-3]-lat[-2]
+        dlon = lon[2]-lon[1]
+        grid_mapping.GeoTransform = f"{lon[0] - dlon/2.} {dlon} 0 {lat[-1] + dlat/2.} 0 {dlat} "
+        print(grid_mapping.GeoTransform)
+
+        spatial_ref = osr.SpatialReference()
+        spatial_ref.ImportFromEPSG(4326)
+        wkt = spatial_ref.ExportToWkt()
+        grid_mapping.spatial_ref = wkt
+
+    output_dataset.sync()
+    output_dataset.close()
+
+
+
+
+
+
+
+if __name__ == "__main__":
+    main()
+