Skip to content

Commit

Permalink
Merge pull request #5 from emit-sds/develop
Browse files Browse the repository at this point in the history
Merge develop into main for v1.2.0
  • Loading branch information
winstonolson authored Jul 3, 2024
2 parents 67abd78 + 4c9eaf4 commit 2e1ad6f
Show file tree
Hide file tree
Showing 6 changed files with 468 additions and 26 deletions.
12 changes: 12 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,22 @@ 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.2.0](https://github.com/emit-sds/emit-sds-l3/compare/v1.1.0...v1.2.0)

> 3 July 2024
- 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)
- add simple vegetation adjustment script for aggregation [`5a26402`](https://github.com/emit-sds/emit-sds-l3/commit/5a264025dd9dfb55534834decfc83a23449c2941)
- add serial apply_glt [`fbb3458`](https://github.com/emit-sds/emit-sds-l3/commit/fbb34582d5e405b77de4e49a8b1f5eb6ee44b49f)
- add support for counting number of revisits [`13d0eed`](https://github.com/emit-sds/emit-sds-l3/commit/13d0eeda2c9dbad5261c4092e183647cd3457409)

#### [v1.1.0](https://github.com/emit-sds/emit-sds-l3/compare/v1.0.0...v1.1.0)

> 6 June 2022
- Merge develop to main for v1.1.0 [`#1`](https://github.com/emit-sds/emit-sds-l3/pull/1)
- removing unmixing.....see emit-sds/SpectralUnmixing [`9d5fa7c`](https://github.com/emit-sds/emit-sds-l3/commit/9d5fa7c48ea669f0eeab3c95fc107d61a5d6f313)
- updated aggregator demo for dynamic load, uncertainty calcs, and netcdf conversion [`cc5241c`](https://github.com/emit-sds/emit-sds-l3/commit/cc5241c4fbf89bba6b8cd75a20d87233322a5e66)
- add license [`247f9b3`](https://github.com/emit-sds/emit-sds-l3/commit/247f9b32c4490ea40cb4c99f30f6569985310dde)
Expand Down
19 changes: 19 additions & 0 deletions CITATION.cff
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
cff-version: 1.2.0
title: emit-sds-l3
message: >-
This codebase provides an interface to executing L3 EMIT
SDS code.
type: software
authors:
- given-names: Philip G Brodrick
affiliation: >-
Jet Propulsion Laboratory, California Institute of
Technology
orcid: 'https://orcid.org/0000-0001-9497-7661'
- given-names: Winston Olson-Duvall
affiliation: >-
Jet Propulsion Laboratory, California Institute of
Technology
orcid: 'https://orcid.org/0000-0002-4210-0283'

license: Apache-2.0
6 changes: 4 additions & 2 deletions apply_glt.py
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,9 @@ def main():


if args.mosaic:
rawspace_files = np.squeeze(np.array(pd.read_csv(args.rawspace_file, header=None)))
rawspace_files = open(args.rawspace_file,'r').readlines()
rawspace_files = [x.strip() for x in rawspace_files]

# TODO: make this check more elegant, should run, catch all files present exception, and proceed
if args.run_with_missing_files is False:
emit_utils.file_checks.check_raster_files(rawspace_files, map_space=False)
Expand Down Expand Up @@ -211,7 +213,7 @@ def apply_mosaic_glt_line(glt_filename: str, output_filename: str, rawspace_file
#glt_line = glt_dataset.ReadAsArray(0, line_index, glt_dataset.RasterXSize, 1)
#glt_line = glt[0][:,line_index:line_index+1, :]

glt_line = np.squeeze(glt[line_index,...]).copy().astype(int)
glt_line = np.squeeze(glt[line_index,...]).copy().astype(int)[...,:3]
valid_glt = np.all(glt_line != GLT_NODATA_VALUE, axis=-1)

glt_line[valid_glt,1] = np.abs(glt_line[valid_glt,1])
Expand Down
141 changes: 141 additions & 0 deletions apply_glt_serial.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,141 @@
"""
Apply a (possibly multi-file) per-pixel spatial reference, in serial (rayless).
Author: Philip G. Brodrick, [email protected]
"""


import argparse
import numpy as np
import pandas as pd
import os
from osgeo import gdal
from spectral.io import envi
import emit_utils.file_checks

from emit_utils.file_checks import envi_header

def _write_bil_chunk(dat, outfile, line, shape, dtype = 'float32'):
"""
Write a chunk of data to a binary, BIL formatted data cube.
Args:
dat: data to write
outfile: output file to write to
line: line of the output file to write to
shape: shape of the output file
dtype: output data type
Returns:
None
"""
outfile = open(outfile, 'rb+')
outfile.seek(line * shape[1] * shape[2] * np.dtype(dtype).itemsize)
outfile.write(dat.astype(dtype).tobytes())
outfile.close()



def single_image_ortho(img_dat, glt, img_ind=None, glt_nodata_value=0):
"""Orthorectify a single image
Args:
img_dat (array like): raw input image
glt (array like): glt - 2 band 1-based indexing for output file(x, y)
img_ind (int): index of image in glt (if mosaic - otherwise ignored)
glt_nodata_value (int, optional): Value from glt to ignore. Defaults to 0.
Returns:
array like: orthorectified version of img_dat
"""
outdat = np.zeros((glt.shape[0], glt.shape[1], img_dat.shape[-1])) - 9999
valid_glt = np.all(glt != glt_nodata_value, axis=-1)

# Only grab data from the correct image, if this is a mosaic
if glt.shape[2] >= 3:
valid_glt[glt[:,:,2] != img_ind] = False

if np.sum(valid_glt) == 0:
return outdat, valid_glt

glt[valid_glt] -= 1 # account for 1-based indexing
outdat[valid_glt, :] = img_dat[glt[valid_glt, 1], glt[valid_glt, 0], :]
return outdat, valid_glt


def main(input_args=None):
parser = argparse.ArgumentParser(description="Robust MF")
parser.add_argument('glt_file', type=str, metavar='GLT', help='path to glt image')
parser.add_argument('raw_file', type=str, metavar='RAW', help='path to raw image')
parser.add_argument('out_file', type=str, metavar='OUTPUT', help='path to output image')
parser.add_argument('--mosaic', action='store_true')
parser.add_argument('--glt_nodata', type=float, default=0)
parser.add_argument('--run_with_missing_files', action='store_true')
parser.add_argument('-b', type=int, nargs='+',default=[-1])
args = parser.parse_args(input_args)


glt_dataset = envi.open(envi_header(args.glt_file))
glt = glt_dataset.open_memmap(writeable=False, interleave='bip').copy().astype(int)
del glt_dataset
glt_dataset = gdal.Open(args.glt_file)

if args.mosaic:
rawspace_files = [x.strip() for x in open(args.raw_file).readlines()]
# TODO: make this check more elegant, should run, catch all files present exception, and proceed
if args.run_with_missing_files is False:
emit_utils.file_checks.check_raster_files(rawspace_files, map_space=False)
# TODO: check that all rawspace files have same number of bands
else:
emit_utils.file_checks.check_raster_files([args.raw_file], map_space=False)
rawspace_files = [args.raw_file]

ort_img = None
band_names = None
for _rf, rawfile in enumerate(rawspace_files):
print(f'{_rf+1}/{len(rawspace_files)}')
if os.path.isfile(envi_header(rawfile)) and os.path.isfile(rawfile):

# Don't load image data unless we have to
if args.mosaic:
if np.sum(glt[:,:,2] == _rf+1) == 0:
continue

img_ds = envi.open(envi_header(rawfile))
inds = None
if args.b[0] == -1:
inds = np.arange(int(img_ds.metadata['bands']))
else:
inds = np.array(args.b)
img_dat = img_ds.open_memmap(writeable=False, interleave='bip')[...,inds].copy()

if band_names is None and 'band names' in envi.open(envi_header(rawfile)).metadata.keys():
band_names = np.array(envi.open(envi_header(rawfile)).metadata['band names'],dtype=str)[inds].tolist()

if ort_img is None:
ort_img, _ = single_image_ortho(img_dat, glt, img_ind=_rf+1, glt_nodata_value=args.glt_nodata)
else:
ort_img_update, valid_glt = single_image_ortho(img_dat, glt, img_ind=_rf+1, glt_nodata_value=args.glt_nodata)
ort_img[valid_glt, :] = ort_img_update[valid_glt, :]

# Build output dataset
driver = gdal.GetDriverByName('ENVI')
driver.Register()

#TODO: careful about output datatypes / format
outDataset = driver.Create(args.out_file, glt.shape[1], glt.shape[0],
ort_img.shape[-1], gdal.GDT_Float32, options=['INTERLEAVE=BIL'])
outDataset.SetProjection(glt_dataset.GetProjection())
outDataset.SetGeoTransform(glt_dataset.GetGeoTransform())
for _b in range(1, ort_img.shape[-1]+1):
outDataset.GetRasterBand(_b).SetNoDataValue(-9999)
if band_names is not None:
outDataset.GetRasterBand(_b).SetDescription(band_names[_b-1])
del outDataset

_write_bil_chunk(ort_img.transpose((0,2,1)), args.out_file, 0, (glt.shape[0], ort_img.shape[-1], glt.shape[1]))




if __name__ == '__main__':
main()


101 changes: 77 additions & 24 deletions build_mosaic_glt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,11 @@ function main()
add_argument!(parser, "--criteria_mode", type = String, default = "distance", help = "Band-ordering criteria mode. Options are min or max (require criteria file), or distance (uses closest point)")
add_argument!(parser, "--criteria_band", type = Int64, default = 1, help = "band of criteria file to use")
add_argument!(parser, "--criteria_file_list", type = String, help = "file(s) to be used for criteria")
add_argument!(parser, "--mask_file_list", type = String, default = nothing, help = "file(s) to be used for mask")
add_argument!(parser, "--mask_band", type = Int64, default = 8, help = "band of mask file to use")
add_argument!(parser, "--target_extent_ul_lr", type = Float64, nargs=4, help = "extent to build the mosaic of")
add_argument!(parser, "--mosaic", type = Int32, default=1, help = "treat as a mosaic")
add_argument!(parser, "--maxlen_file", type = String, default=nothing, help = "max length reference file")
add_argument!(parser, "--output_epsg", type = Int32, default=4326, help = "epsg to write to destination")
add_argument!(parser, "--log_file", type = String, default = nothing, help = "log file to write to")
args = parse_args(parser)
Expand Down Expand Up @@ -54,6 +57,22 @@ function main()
# TODO: add check to make sure criteria file dimensions match igm file dimensions
end

if !isnothing(args.mask_file_list)
if args.mosaic == 1
mask_files = readdlm(args.mask_file_list, String)
else
mask_files = [args.mask_file_list]
end
end

if !isnothing(args.maxlen_file)
if args.mosaic == 1
maxlen_files = readdlm(args.maxlen_file, String)
else
maxlen_files = [args.maxlen_file]
end
end

if length(args.target_extent_ul_lr) > 0
ullr = args.target_extent_ul_lr
min_x = ullr[1]
Expand All @@ -78,7 +97,7 @@ function main()

@info "Output Image Size (x,y): $x_size_px, $y_size_px. Creating output dataset."
if args.mosaic == 1
output_bands = 3
output_bands = 4
else
output_bands = 2
end
Expand All @@ -93,34 +112,56 @@ function main()
grid[..,2] = fill(1,y_size_px,x_size_px) .* LinRange(max_y + args.target_resolution[2]/2,max_y + args.target_resolution[2] * (1/2 + y_size_px - 1), y_size_px)[:,[CartesianIndex()]]

@info "Create GLT."
best = fill(1e12, y_size_px, x_size_px, 4)
best = fill(1e12, y_size_px, x_size_px, 5)
if args.criteria_mode == "max"
best = best .* -1
end
best[..,1:3] .= -9999
best[..,5] .= 0

max_offset_distance = sqrt(sum(args.target_resolution.^2))*3
pixel_buffer_window = 1

total_found = 0
lk = ReentrantLock()
for (file_idx, igm_file) in enumerate(igm_files)
@info "$igm_file"
dataset = ArchGDAL.read(igm_file)
dataset = ArchGDAL.read(igm_file,alloweddrivers =["ENVI"])
igm = PermutedDimsArray(ArchGDAL.read(dataset), (2,1,3))
if minimum(igm[..,1]) > grid[1,end-1,1] || maximum(igm[..,1]) < grid[1,1,1] ||
minimum(igm[..,2]) > grid[1,1,2] || maximum(igm[..,2]) < grid[end-1,1,2]
#println(minimum(igm[..,1]), " > ", grid[1,end-1,1], " ", maximum(igm[..,1]), " < ", grid[1,1,1])
#println(minimum(igm[..,2]), " > ", grid[1,1,2], " ", maximum(igm[..,2]), " < ", grid[end-1,1,2])
continue
else
println("Entering")
end

#if minimum(igm[..,1]) > grid[1,end-1,1] || maximum(igm[..,1]) < grid[1,1,1] ||
# minimum(igm[..,2]) > grid[1,1,2] || maximum(igm[..,2]) < grid[end-1,1,2]
# #println(minimum(igm[..,1]), " > ", grid[1,end-1,1], " ", maximum(igm[..,1]), " < ", grid[1,1,1])
# #println(minimum(igm[..,2]), " > ", grid[1,1,2], " ", maximum(igm[..,2]), " < ", grid[end-1,1,2])
# continue
#else
# println("Entering")
#end
if args.criteria_mode != "distance"
criteria_dataset = ArchGDAL.read(criteria_files[file_idx])
cffi = criteria_files[file_idx]
@debug "$cffi"
criteria_dataset = ArchGDAL.read(criteria_files[file_idx],alloweddrivers =["ENVI"])
criteria = PermutedDimsArray(ArchGDAL.read(criteria_dataset, args.criteria_band), (2,1))
end
for _y=1:size(igm)[1]
for _x=1:size(igm)[2]
if !isnothing(args.mask_file_list)
mffi = mask_files[file_idx]
@debug "$mffi"
mask_dataset = ArchGDAL.read(mask_files[file_idx],alloweddrivers =["ENVI"])
mask = PermutedDimsArray(ArchGDAL.read(mask_dataset, args.mask_band), (2,1))
end
if !isnothing(args.maxlen_file)
maxlen_ds = ArchGDAL.readraster(maxlen_files[file_idx],alloweddrivers=["ENVI"])
maxlines = size(maxlen_ds)[2]
else
maxlines = size(igm)[1]
end

Threads.@threads for _y=1:min(size(igm)[1],maxlines)
Threads.@threads for _x=1:size(igm)[2]
if !isnothing(args.mask_file_list)
if mask[_y, _x] > 0
continue
end
end
pt = igm[_y,_x,1:2]
closest_t = Array{Int64}([round((pt[2] - grid[1,1,2]) / args.target_resolution[2]),
round((pt[1] - grid[1,1,1]) / args.target_resolution[1]) ]) .+ 1
Expand All @@ -131,31 +172,43 @@ function main()
closest[1] = closest_t[1] + xbuffer
closest[2] = closest_t[2] + ybuffer


if closest[1] < 1 || closest[2] < 1 || closest[1] > size(grid)[1] || closest[2] > size(grid)[2]
continue
end
dist = sum((grid[closest[1],closest[2],:] - pt).^2)

if dist < max_offset_distance

best[closest[1], closest[2], 5] += 1
if args.criteria_mode in ["distance", "min"]
if args.criteria_mode == "distance"
current_crit = dist
else
current_crit = criteria[_y, _x]
end

if current_crit < best[closest[1], closest[2], 4]
best[closest[1], closest[2], 1:3] = [_x, _y, file_idx]
best[closest[1], closest[2], 4] = current_crit
lock(lk)
try
if current_crit < best[closest[1], closest[2], 4]
best[closest[1], closest[2], 1:3] = [_x, _y, file_idx]
best[closest[1], closest[2], 4] = current_crit
end
finally
unlock(lk)
end

elseif args.criteria_mode == "max"
current_crit = criteria[_y, _x]
if current_crit > best[closest[1], closest[2], 4]
best[closest[1], closest[2], 1:3] = [_x, _y, file_idx]
best[closest[1], closest[2], 4] = current_crit
lock(lk)
try
if current_crit > best[closest[1], closest[2], 4]
best[closest[1], closest[2], 1:3] = [_x, _y, file_idx]
best[closest[1], closest[2], 4] = current_crit
end
finally
unlock(lk)
end

end
end
end
Expand All @@ -165,9 +218,9 @@ function main()
end
end

println(total_found, " ", sum(best[..,1] .!= -9999), " ", size(best)[1]*size(best)[2])
println(sum(best[..,1] .!= -9999), " ", size(best)[1]*size(best)[2])
if args.mosaic == 1
output = Array{Int32}(permutedims(best[..,1:3], (2,1,3)))
output = Array{Int32}(permutedims(best[..,[1,2,3,5]], (2,1,3)))
else
output = Array{Int32}(permutedims(best[..,1:2], (2,1,3)))
end
Expand Down
Loading

0 comments on commit 2e1ad6f

Please sign in to comment.