diff --git a/demcoreg/dem_mask.py b/demcoreg/dem_mask.py index 7968434..734fa7d 100755 --- a/demcoreg/dem_mask.py +++ b/demcoreg/dem_mask.py @@ -31,11 +31,11 @@ datadir = iolib.get_datadir() -def get_nlcd_fn(yr=2016): +def get_nlcd_fn(yr=2019): """Calls external shell script `get_nlcd.sh` to fetch: Land Use Land Cover (nlcd) grids, 30 m - 2011, 2013 or 2016 (default) + 2011, 2013, 2016, or 2019 (default) http://www.mrlc.gov/nlcd11_leg.php """ @@ -69,6 +69,15 @@ def get_bareground_fn(): #subprocess.call(cmd) return bg_fn +def get_canopyheight_fn(region='NAM'): + canopyheight_fn = os.path.join(datadir, f'canopyheight/Forest_height_2019_{region}.tif') + print(canopyheight_fn) + if not os.path.exists(canopyheight_fn): + cmd = ['get_canopyheight.sh', region] + #subprocess.call(cmd) + sys.exit("Missing canopy height data source. If already downloaded, specify correct datadir. If not, run `%s` to download" % cmd[0]) + return canopyheight_fn + #Download latest global RGI glacier db def get_glacier_poly(): """Calls external shell script `get_rgi.sh` to fetch: @@ -133,6 +142,9 @@ def get_nlcd_mask(nlcd_ds, filter='not_forest', out_fn=None): mask = ~(np.logical_or(np.logical_or((l==41),(l==42)),(l==43))) elif filter == 'not_forest+not_water': mask = ~(np.logical_or(np.logical_or(np.logical_or((l==41),(l==42)),(l==43)),(l==11))) + elif filter == 'barren+flatlowveg': + # rock+ice+shrub+scrub+grass+sedge+pasture+wetlands was too long + mask = np.logical_or(np.logical_or(np.logical_or(np.logical_or(np.logical_or(np.logical_or(np.logical_or(np.logical_or((l==12),(l==31)),(l==51)),(l==52)),(l==71)),(l==72)),(l==81)),(l==90)),(l==95)) else: print("Invalid mask type") mask = None @@ -160,6 +172,36 @@ def get_bareground_mask(bareground_ds, bareground_thresh=60, out_fn=None): l = None return mask +# Create canopy height mask +def get_canopyheight_mask(canopyheight_ds, filter='ground', out_fn=None): + """Generate raster mask for specified global forest canopy height vegetation class filter + """ + print("Loading global forest canopy height") + b = canopyheight_ds.GetRasterBand(1) + l = b.ReadAsArray() + print("Filtering canopy height with: %s" % filter) + # Original canopy height products use 103 as ndv + #0–60 - canopy height [m] + #101 - water + #102 - snow/ice + if filter == 'ground': + mask = (l==1) + elif filter == 'water': + mask = (l==101) + elif filter == 'snow_ice': + mask = (l==102) + elif filter == 'water+snow_ice': + mask = np.logical_or((l==101),(l==102)) + else: + print("Invalid mask type") + mask = None + #Write out original data + if out_fn is not None: + print("Writing out %s" % out_fn) + iolib.writeGTiff(l, out_fn, canopyheight_ds) + l = None + return mask + def get_snodas_ds(dem_dt, code=1036): """Function to fetch and process SNODAS snow depth products for input datetime @@ -269,6 +311,11 @@ def get_modscag_fn_list(dem_dt, tile_list=('h08v04', 'h09v04', 'h10v04', 'h08v05 import requests from bs4 import BeautifulSoup auth = iolib.get_auth() + # from requests.auth import HTTPDigestAuth + # uname='' + # pw='' + # auth = HTTPDigestAuth(uname, pw) + pad_days = timedelta(days=pad_days) dt_list = timelib.dt_range(dem_dt-pad_days, dem_dt+pad_days+timedelta(1), timedelta(1)) @@ -282,9 +329,16 @@ def get_modscag_fn_list(dem_dt, tile_list=('h08v04', 'h09v04', 'h10v04', 'h08v05 #If we already have a vrt and it contains all of the necessary tiles if os.path.exists(out_vrt_fn): vrt_ds = gdal.Open(out_vrt_fn) - if np.all([np.any([tile in sub_fn for sub_fn in vrt_ds.GetFileList()]) for tile in tile_list]): + # hack to avoid tile checking on M1 + checktiles=False + if checktiles: + if np.all([np.any([tile in sub_fn for sub_fn in vrt_ds.GetFileList()]) for tile in tile_list]): + out_vrt_fn_list.append(out_vrt_fn) + continue + else: out_vrt_fn_list.append(out_vrt_fn) continue + #Otherwise, download missing tiles and rebuild #Try to use historic products modscag_fn_list = [] @@ -298,7 +352,7 @@ def get_modscag_fn_list(dem_dt, tile_list=('h08v04', 'h09v04', 'h10v04', 'h08v05 modscag_url_fn = [] if r.ok: parsed_html = BeautifulSoup(r.content, "html.parser") - modscag_url_fn = parsed_html.findAll(text=re.compile('%s.*snow_fraction.tif' % tile)) + modscag_url_fn = parsed_html.find_all(href=re.compile(f'{tile}.*?snow_fraction.tif')) if not modscag_url_fn: #Couldn't find historic, try to use real-time products modscag_url_str = 'https://snow-data.jpl.nasa.gov/modscag/%Y/%j/' @@ -307,7 +361,7 @@ def get_modscag_fn_list(dem_dt, tile_list=('h08v04', 'h09v04', 'h10v04', 'h08v05 r = requests.get(modscag_url_base, auth=auth) if r.ok: parsed_html = BeautifulSoup(r.content, "html.parser") - modscag_url_fn = parsed_html.findAll(text=re.compile('%s.*snow_fraction.tif' % tile)) + modscag_url_fn = parsed_html.find_all(href=re.compile(f'{tile}.*?snow_fraction.tif')) if not modscag_url_fn: print("Unable to fetch MODSCAG for %s" % dt) else: @@ -315,12 +369,14 @@ def get_modscag_fn_list(dem_dt, tile_list=('h08v04', 'h09v04', 'h10v04', 'h08v05 #Now extract actual tif filenames to fetch from html parsed_html = BeautifulSoup(r.content, "html.parser") #Fetch all tiles - modscag_url_fn = parsed_html.findAll(text=re.compile('%s.*snow_fraction.tif' % tile)) - if modscag_url_fn: + modscag_url_fn = parsed_html.find_all(href=re.compile(f'{tile}.*?snow_fraction.tif')) + if modscag_url_fn: modscag_url_fn = modscag_url_fn[0] + modscag_url_fn=modscag_url_fn.get('href') modscag_url = os.path.join(modscag_url_base, modscag_url_fn) print(modscag_url) - modscag_fn = os.path.join(outdir, os.path.split(modscag_url_fn)[-1]) + #modscag_fn = os.path.join(outdir, os.path.split(modscag_url_fn)[-1]) + modscag_fn = os.path.join(outdir, modscag_url_fn) if not os.path.exists(modscag_fn): iolib.getfile2(modscag_url, auth=auth, outdir=outdir) modscag_fn_list.append(modscag_fn) @@ -432,6 +488,8 @@ def get_mask(dem_ds, mask_list, dem_fn=None, writeout=False, outdir=None, args=N if dem_fn is not None: #Extract DEM timestamp dem_dt = timelib.fn_getdatetime(dem_fn) + if dem_dt is None and args.dem_dt is not None: + dem_dt=datetime.strptime(args.dem_dt, "%Y%m%d") out_fn_base = os.path.join(outdir, os.path.splitext(os.path.split(dem_fn)[-1])[0]) if args is None: @@ -504,17 +562,18 @@ def get_mask(dem_ds, mask_list, dem_fn=None, writeout=False, outdir=None, args=N print("SNODAS grid for input location and timestamp is empty") #These tiles cover CONUS - #tile_list=('h08v04', 'h09v04', 'h10v04', 'h08v05', 'h09v05') + tile_list=('h08v04', 'h09v04', 'h10v04', 'h08v05', 'h09v05') if 'modscag' in mask_list and args.modscag_thresh > 0: modscag_min_dt = datetime(2000,2,24) if dem_dt < modscag_min_dt: print("Warning: DEM timestamp (%s) is before earliest MODSCAG timestamp (%s)" \ % (dem_dt, modscag_min_dt)) else: - tile_list = get_modis_tile_list(dem_ds) + # tile_list = get_modis_tile_list(dem_ds) print(tile_list) pad_days=7 modscag_fn_list = get_modscag_fn_list(dem_dt, tile_list=tile_list, pad_days=pad_days) + print('modscag_fn_list is', modscag_fn_list) if modscag_fn_list: modscag_ds = proc_modscag(modscag_fn_list, extent=dem_ds, t_srs=dem_ds) modscag_ds_warp = warplib.memwarp_multi([modscag_ds,], res=dem_ds, extent=dem_ds, t_srs=dem_ds, r='cubicspline')[0] @@ -528,7 +587,7 @@ def get_mask(dem_ds, mask_list, dem_fn=None, writeout=False, outdir=None, args=N modscag_mask = (modscag_fsca.filled(0) >= args.modscag_thresh) modscag_mask = ~(modscag_mask) if writeout: - out_fn = os.path.splitext(out_fn)[0]+'_mask.tif' + out_fn = os.path.splitext(out_fn)[0]+f'_{args.modscag_thresh}thresh_mask.tif' print("Writing out %s" % out_fn) iolib.writeGTiff(modscag_mask, out_fn, src_ds=dem_ds) newmask = np.logical_and(modscag_mask, newmask) @@ -545,6 +604,24 @@ def get_mask(dem_ds, mask_list, dem_fn=None, writeout=False, outdir=None, args=N iolib.writeGTiff(toa_mask, out_fn, src_ds=dem_ds) newmask = np.logical_and(toa_mask, newmask) + if 'canopy' in mask_list: + # Will need to think about how to warp this... + canopyheight_ds = gdal.Open(get_canopyheight_fn()) + canopyheight_ds_warp = warplib.memwarp_multi([canopyheight_ds,], res=dem_ds, extent=dem_ds, t_srs=dem_ds, r='cubicspline')[0] + out_fn = None + if writeout: + out_fn = out_fn_base+'_canopy.tif' + canopy_mask = get_canopyheight_mask(canopyheight_ds_warp, filter=args.canopy_filter, out_fn=out_fn) + # ground = iolib.ds_getma(canopy_ds_warp) + # canopy_mask = (ground.filled(0) > 1) + # canopy_mask = ~(canopy_mask) + + if writeout: + out_fn = os.path.splitext(out_fn)[0]+'_mask.tif' + print("Writing out %s" % out_fn) + iolib.writeGTiff(canopy_mask, out_fn, src_ds=dem_ds) + newmask = np.logical_and(canopy_mask, newmask) + if False: #Filter based on expected snowline #Simplest approach uses altitude cutoff @@ -567,7 +644,7 @@ def get_mask(dem_ds, mask_list, dem_fn=None, writeout=False, outdir=None, args=N return newmask #Can add "mask_list" argument, instead of specifying individually -mask_choices = ['toa', 'snodas', 'modscag', 'bareground', 'glaciers', 'nlcd', 'none'] +mask_choices = ['toa', 'snodas', 'modscag', 'bareground', 'glaciers', 'nlcd', 'canopy', 'none'] def getparser(): parser = argparse.ArgumentParser(description="Identify control surfaces for DEM co-registration") parser.add_argument('dem_fn', type=str, help='DEM filename') @@ -584,9 +661,13 @@ def getparser(): parser.add_argument('--bareground_thresh', type=float, default=60, help='Percent bareground threshold (default: %(default)s%%, valid range 0-100), mask greater than this value (only relevant for global bareground data)') parser.add_argument('--glaciers', action='store_true', help="Mask glacier polygons") parser.add_argument('--nlcd', action='store_true', help="Enable NLCD LULC filter (for CONUS)") - nlcd_filter_choices = ['rock', 'rock+ice', 'rock+ice+water', 'not_forest', 'not_forest+not_water', 'none'] + nlcd_filter_choices = ['rock', 'rock+ice', 'rock+ice+water', 'not_forest', 'not_forest+not_water', 'barren+flatlowveg', 'none'] parser.add_argument('--nlcd_filter', type=str, default='not_forest', choices=nlcd_filter_choices, help='Preserve these NLCD pixels (default: %(default)s)') + parser.add_argument('--canopy', action='store_true', help="Enable global canopy height filter (for CONUS)") + canopy_filter_choices = ['ground', 'water', 'snow_ice', 'water+snow_ice', 'none'] + parser.add_argument('--canopy_filter', type=str, default='ground', choices=canopy_filter_choices, help='Preserve these canopy height pixels (default: %(default)s)') parser.add_argument('--dilate', type=int, default=None, help='Dilate mask with this many iterations (default: %(default)s)') + parser.add_argument('--dem_dt', type=str, default=None, help='Manually input DEM datetime as YYYYMMDD') return parser def main(): @@ -600,6 +681,7 @@ def main(): if args.bareground: mask_list.append('bareground') if args.glaciers: mask_list.append('glaciers') if args.nlcd: mask_list.append('nlcd') + if args.canopy: mask_list.append('canopy') if not mask_list: parser.print_help() @@ -638,12 +720,16 @@ def main(): min_validpx_count = 100 min_validpx_std = 10 validpx_count = newdem.count() - validpx_std = newdem.std() + # validpx_std = newdem.std() + validpx_std = np.nanstd(newdem) #.std() print("%i valid pixels in masked output tif to be used as ref" % validpx_count) print("%0.2f std in masked output tif to be used as ref" % validpx_std) #if (validpx_count > min_validpx_count) and (validpx_std > min_validpx_std): if (validpx_count > min_validpx_count): - out_fn = os.path.join(args.outdir, os.path.splitext(os.path.split(dem_fn)[-1])[0]+'_ref.tif') + if args.modscag: + out_fn = os.path.join(args.outdir, os.path.splitext(os.path.split(dem_fn)[-1])[0]+f'_{args.modscag_thresh}thresh_ref.tif') + else: + out_fn = os.path.join(args.outdir, os.path.splitext(os.path.split(dem_fn)[-1])[0]+'_ref.tif') print("Writing out %s" % out_fn) iolib.writeGTiff(newdem, out_fn, src_ds=dem_ds) else: diff --git a/demcoreg/dem_stable_zalign.py b/demcoreg/dem_stable_zalign.py new file mode 100755 index 0000000..6ee805a --- /dev/null +++ b/demcoreg/dem_stable_zalign.py @@ -0,0 +1,504 @@ +#! /usr/bin/env python + +import sys +import os +import argparse +import subprocess + +import numpy as np +import matplotlib.pyplot as plt + +from pygeotools.lib import warplib, iolib, malib, geolib +from imview.lib import pltlib + +from pathlib import PurePath +import rasterio as rio +from rasterio.enums import Resampling + +def load(fn, name, masked=True, squeeze=True, dtype=np.float32, chunks='auto'): + """Function to open files with rasterio as xarray DataArrays + + Parameters + ------------- + thisDir: str + directory path to search + fn_pattern: str + regex pattern to match files + verbose: boolean + print filenames + + Returns + ------------- + fns: list + list of filenames matched and sorted + """ + import rioxarray + if squeeze: + arr=rioxarray.open_rasterio(fn, masked=masked, default_name=name, chunks=chunks).squeeze(dim='band', drop=True) + else: + arr=rioxarray.open_rasterio(fn, masked=masked, default_name=name, chunks=chunks) + + if dtype: arr=arr.astype(dtype) + + return arr + +def coerce_float(a, delimiter='+', verbose=False): + """Docstring short description + Parameters + ------------- + thisDir: str + directory path to search + fn_pattern: str + regex pattern to match files + verbose: boolean + print filenames + + Returns + ------------- + fns: list + list of filenames matched and sorted + """ + if verbose: print(a) + try: + a=float(a) + except: + a = a.split('+')[-1] + if verbose: print(a) + a=float(a) + if verbose: print(a) + return a + +def extract_shifts(align_fn, ref_fn, verbose=False): + """Extract x, y, and z adjustments based on input align fn and coerce to float + + Parameters + ------------- + align_fn: str + aligned filename + ref_fn: str + reference filename + verbose: boolean + print statement checks + + Returns + ------------- + x, y, z : tuple + tuple of float shifts + """ + x = os.path.basename(align_fn).split(f'{os.path.basename(ref_fn)[:-4]}')[-1].split('_x')[1].split('_')[0] + y = os.path.basename(align_fn).split(f'{os.path.basename(ref_fn)[:-4]}')[-1].split('_y')[1].split('_')[0] + z = os.path.basename(align_fn).split(f'{os.path.basename(ref_fn)[:-4]}')[-1].split('_z')[1].split('_')[0] + if verbose: print(x, y, z) + x = coerce_float(x) + y = coerce_float(y) + z = coerce_float(z) + if verbose: print(x, y, z) + if verbose: print(type(x), type(y), type(z)) + return x, y, z + +def zshiftstr(zadjust): + """Pre-pend sign for zshift >=0, add "+" and retain hundredths place + + Parameters + ------------- + zadjust: float + dz to stringify + + Returns + ------------- + zstr: str + string with prepended sign + """ + if zadjust >= 0: zstr=f'+{zadjust:.2f}' + else: zstr = f'{zadjust:.2f}' + print(zstr) + return zstr + +def find_mode(arr, dec=2, verbose=False): + '''By default, omits nans if present and computes mode over entire array. + Can return multiple modes. + If arr dtype is float, will round values to the specified decimals + e.g., 1 = tenths (0.1), 2 = hundredths (0.01) + + Parameters + ------------- + arr: np.array + array from which to derive mode + dec: int + decimal places to round + verbose: boolean + print mode value and count + + Returns + ------------- + modes: list + list of identified modes + ''' + # first remove nans + arr = arr[~np.isnan(arr)] + + # then round the values + arr = arr * 10**dec + arr = arr.astype(int) + + # then compute mode(s) by using a dict and a counter loop approach + # this could take awhile for large arrays + # https://www.nobledesktop.com/learn/python/mode-in-python + vals = {} + for v in arr: + if not v in vals: + vals[v] = 1 + else: + vals[v]+=1 + + # this will return all the modes + modes = [g for g,l in vals.items() if l==max(vals.values())] + if verbose: + _ = [print(f'{m / 10**dec}: {vals[m]}') for m in modes] + modes = [m / 10**dec for m in modes] + return modes + +def getparser(): + parser = argparse.ArgumentParser(description="Perform two-stage co-registration using stable surface z shift", \ + formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument('ref_fn', type=str, help='Reference DEM filename') + parser.add_argument('stable_fn', type=str, help='Stable surface masked referenced DEM filename to be shifted') + parser.add_argument('align_fn', type=str, help='Full surface aligned DEM filename to be stable surface median z shifted') + parser.add_argument('-new_align_fn', default=None, type=str, help='Updated full surface two stage aligned DEM filename to output') + parser.add_argument('-outdir', default=None, help='Output directory') + parser.add_argument('-realrun', default=True, type=bool, help='Run calculations and write out files') + parser.add_argument('-v', default=False, type=bool, help='Increase verbosity') + return parser + +def main(argv=None): + parser = getparser() + args = parser.parse_args() + + #Should check that files exist + ref_fn = args.ref_fn + stable_fn = args.stable_fn + align_fn = args.align_fn + updated_align_fn = args.new_align_fn + realrun = args.realrun + verbose = args.v + outdir = args.outdir + + # consider adding these as built in args + overwrite = False + resampling = Resampling.average + + if outdir is None: + outdir = PurePath(align_fn).parents[0] + elif not os.path.exists(outdir): + os.makedirs(outdir) + + outprefix = f'{PurePath(align_fn).stem.split("_nuth_")[0]}' + + print(f'\nCurrent directory: {os.getcwd()}') + print(f'Running: dem_stable_zalign.py') + + print(f"\nReference: {ref_fn}") + print(f"Stable reference: {stable_fn}") + print(f"Full surface aligned DEM: {align_fn}") + + # Extract shifts from align_fn based on input reference + # Note the horizontal shifts (xy) from full shift alignment png + x, y, z = extract_shifts(align_fn, ref_fn=ref_fn, verbose=verbose) + # print(x, y, z) + + # Load the aligned and the stable surface tif + align_arr = load(align_fn, PurePath(align_fn).stem) + stable_arr = load(stable_fn, PurePath(stable_fn).stem) + + # Explicitly write out no data values + align_arr.rio.write_nodata(np.nan, inplace=True) + stable_arr.rio.write_nodata(np.nan, inplace=True) + + if verbose: + print('\nInput array resolutions:') + print(f'align: {align_arr.rio.resolution()[0]}') + print(f'stable ref: {stable_arr.rio.resolution()[0]}') + + # resample the stable surface reference to the aligned array + # matching extents and the input aligned resolution + stable_arr_reproj = stable_arr.rio.reproject_match(align_arr, + resampling=resampling, + nodata=np.nan) + if verbose: print(f'stable ref reproject-matched: {stable_arr_reproj.rio.resolution()[0]}') + + # if verbose: print(f'Reproject matched check: {stable_arr_reproj.shape}, {align_arr.shape}') + + # create binary mask to remove smooshed non-stable surfaces + import copy + stable_binary = copy.deepcopy(stable_arr) + + stable_binaryval = stable_binary.values + if verbose: print(f'Starting NaNcount: {np.isnan(stable_binaryval).sum()}') + stable_binaryval[~np.isnan(stable_binaryval)] = 1 + stable_binaryval[np.isnan(stable_binaryval)] = 0 + stable_binary.values = stable_binaryval + if verbose: print(f'Ending NaNcount after binarization: {np.isnan(stable_binary.values).sum()}') + + # Reproject the binary stable arr to match align_arr + stable_binary_arr_reproj = stable_binary.rio.reproject_match(align_arr, + resampling=resampling, + nodata=np.nan) + if verbose: print(f'Shapes check: {stable_binary_arr_reproj.shape}, {stable_arr_reproj.shape}, {align_arr.shape}') + + # Set pixel values that do not meet threshold to nan + strict_thresh = 1 + if verbose: print(f'Starting NaNcount before thresholding: {np.isnan(stable_arr_reproj.values).sum()}') + stable_arr_reprojval = stable_arr_reproj.values + stable_binary_arr_reprojval = stable_binary_arr_reproj.values + stable_arr_reprojval[stable_binary_arr_reprojval < strict_thresh] = np.nan + stable_arr_reproj.values = stable_arr_reprojval + if verbose: print(f'Ending NaNcount before thresholding: {np.isnan(stable_arr_reproj.values).sum()}') + + # Undo the z shift from the warped full surface align.tif by adding the + # opposite zshift to both the horizontal and the full array + if verbose: print('\nRe-setting z shifts') + align_warp_xyalone = align_arr + -z + align_arr = align_arr + -z + + # Compute the median difference over stable surface areas using the + # matched stable surface and the undone horizontally aligned array + # this is done to match the align.tif (coarser) array + print('\nComputing median difference over stable surfaces') + stable_diff = align_warp_xyalone - stable_arr_reproj + # if verbose: print(stable_diff.shape) + + # Calculate the median difference over stable surfaces to be used for z shift + zshift = np.nanmedian(stable_diff.values) + + print(f'The median diff of (src - ref) is: {zshift:.2f} m') + + calcmodes=True + if calcmodes: + # Include this bit to check the diff between med and mode + modezshift = find_mode(stable_diff.values) + if len(modezshift) > 1: + print(f'The mode diffs of (src - ref) is:') + _ = [print(f'{m:.2f} m') for m in modezshift] + else: + print(f'The mode diff of (src - ref) is: {modezshift[0]:.2f} m') + + # sys.exit(1) + # if this is still nan after all the rigamarole, then exit out on error + if np.isnan(zshift): + print('Z shift is nan, check for overlapping stable surfaces') + print("\nReference: %s" % ref_fn) + print("Stable reference: %s" % stable_fn) + print("Full surface aligned DEM: %s" % align_fn) + sys.exit(1) + + # Note the z shift here for writing out file + zstr = f'{-zshift:.2f}' + print(f'z shift string for source adjustment is: {zshiftstr(float(zstr))}') + + # Apply z shift to the xy-aligned array that moves the median difference to 0 + print('\nApplying z shift') + align_warp_fullshift = align_warp_xyalone - zshift + + # Apply the z shift to the full aligned array, not the warped horizontal array + align_arr = align_arr - zshift + + # This diff is also src-ref now + stable_diff_updated = align_warp_fullshift - stable_arr_reproj + + # check that the median diff is now 0 over stable surfaces + print(f'New median diff ≈ 0? \ + {np.isclose(np.nanmedian(stable_diff_updated.values), 0, atol=1e-3)}') + + # Specify filename - this assumes nuth and kaab alignment, could leverage + # the x bit + if updated_align_fn is None: + updated_align_fn = f'{align_fn.split("_nuth_")[0]}_nuth_twostageshift_x{zshiftstr(x)}_y{zshiftstr(y)}_z{zshiftstr(float(zstr))}_align.tif' + if verbose: print(f'\nNew align.tif: {updated_align_fn}') + + # Update attributes and the name of array + align_arr.name = PurePath(updated_align_fn).stem + align_arr.attrs = {'stable_fn': stable_fn} + + # and write it out + if not os.path.exists(updated_align_fn): + print("\nDNE, writing out now...") + print(updated_align_fn) + if realrun: align_arr.rio.to_raster(updated_align_fn) + else: + print('Exists') + + # Specify diff fn + updated_diff_fn = f'{updated_align_fn[:-4]}_diff.tif' + if not os.path.exists(updated_diff_fn): + + # Pull in the full ref_fn + ref_arr = load(ref_fn, PurePath(ref_fn).stem) + ref_arr.rio.write_nodata(np.nan, inplace=True) + + # reproject and match the original align arr + ref_arr_reproj = ref_arr.rio.reproject_match(align_arr, + resampling=resampling, + nodata=np.nan) + if verbose: print(ref_arr_reproj.rio.resolution()) + + # Calculated the updated diff for all surfaces using the updated align_arr + diff_updated = align_arr - ref_arr_reproj + + # and write it out + print("\nDNE, writing out now...") + print(updated_diff_fn) + if realrun: diff_updated.rio.to_raster(updated_diff_fn) + del diff_updated + else: + print('Exists') + + del x, y, z, align_arr, stable_arr, stable_arr_reproj, align_warp_xyalone, \ + stable_diff, zshift, zstr, align_warp_fullshift, stable_diff_updated + + print('\nMoving to stats and plotting section') + align_stats_fn = f'{updated_align_fn[:-4]}_stats.json' + # print(updated_align_fn) + if os.path.exists(align_stats_fn) and not overwrite: + print('Stats file exists, overwrite set to False') + print(align_stats_fn) + else: + # consider adding json stats bits and plotting code here + orig_diff_fn = f'{align_fn.split("_nuth")[0]}_orig_diff.tif' + ref_base = PurePath(ref_fn).stem + # dem_fn = f'{PurePath(align_fn).parents[1]}/{PurePath(align_fn.split("_" + ref_base)[0]).stem}.tif' # worked for stacks + # dem_fn = f'{PurePath(align_fn).parents[0].stem}/{PurePath(align_fn.split("_" + ref_base)[0]).stem}.tif' # for the stacks_clean set up, Pléiades + dem_fn = f'{align_fn.split("_dem_align")[0]}.tif' + print(dem_fn, ref_fn, orig_diff_fn) + src_fn_list = [dem_fn, ref_fn, orig_diff_fn] + + elev_dslist = warplib.memwarp_multi_fn(src_fn_list=src_fn_list, res='mean') + src_dem_clip_ds = elev_dslist[0] + ref_dem_clip_ds = elev_dslist[1] + ref_dem_orig = iolib.ds_getma(ref_dem_clip_ds) + src_dem_orig = iolib.ds_getma(src_dem_clip_ds) + + # add bits for stable surface shifting plot + _, stable_ds = warplib.memwarp_multi_fn([orig_diff_fn, stable_fn], res='first') + stable_orig = iolib.ds_getma(stable_ds) + + #Needed for plotting + ref_dem_hs = geolib.gdaldem_mem_ds(ref_dem_clip_ds, processing='hillshade', returnma=True, computeEdges=True) + src_dem_hs = geolib.gdaldem_mem_ds(src_dem_clip_ds, processing='hillshade', returnma=True, computeEdges=True) + diff_orig = src_dem_orig - ref_dem_orig + diff_orig_stats = malib.get_stats_dict(diff_orig, full=True) + # adapt for stable surface stat plotting only + + align_dslist = warplib.memwarp_multi_fn(src_fn_list=[updated_align_fn, ref_fn], res='mean', dst_ndv=-9999) + src_dem_align = iolib.ds_getma(align_dslist[0]) + ref_dem_align = iolib.ds_getma(align_dslist[1]) + diff_align = src_dem_align - ref_dem_align + diff_align_stats = malib.get_stats_dict(diff_align[~np.isnan(diff_align)], full=True) + + res = geolib.get_res(ref_dem_clip_ds)[0] + + + dx, dy, dz = extract_shifts(updated_align_fn, ref_fn=ref_fn) + dm_total = np.sqrt(dx**2 + dy**2 + dz**2) + + # write out stats jsons + align_stats_fn = f'{updated_align_fn[:-4]}_stats.json' + print(f'Pulling stats for {align_stats_fn}') + align_stats = {} + align_stats['src_fn'] = dem_fn + align_stats['ref_fn'] = ref_fn + align_stats['stable_fn'] = stable_fn + align_stats['updated_align_fn'] = updated_align_fn + align_stats['shift'] = {'dx':dx, 'dy':dy, 'dz':dz, 'dm':dm_total} + align_stats['before'] = diff_orig_stats + align_stats['after'] = diff_align_stats + + import json + with open(align_stats_fn, 'w') as f: + print('Dumping stats json') + if realrun: json.dump(align_stats, f) + + if realrun: + fig_fn = f'{updated_align_fn[:-4]}.png' + + if os.path.exists(fig_fn) and not overwrite: + print('Updated figure png exists') + print(fig_fn) + else: + + f, axa = plt.subplots(2, 4, figsize=(16, 8)) + for ax in axa.ravel()[:-1]: + ax.set_facecolor('k') + pltlib.hide_ticks(ax) + dem_clim = malib.calcperc(ref_dem_orig, (2,98)) + # cmap = plt.cm.jet + cmap = 'cpt_rainbow' # plt.cm.rainbow + + kwargs = {'interpolation':'none'} + axa[0,0].imshow(ref_dem_hs, cmap='gray', **kwargs) + im = axa[0,0].imshow(ref_dem_orig, cmap=cmap, clim=dem_clim, alpha=0.6, **kwargs) + pltlib.add_cbar(axa[0,0], im, arr=ref_dem_orig, clim=dem_clim, label=None) + pltlib.add_scalebar(axa[0,0], res=res) + axa[0,0].set_title('Reference DEM') + + axa[0,1].imshow(src_dem_hs, cmap='gray', **kwargs) + im = axa[0,1].imshow(src_dem_orig, cmap=cmap, clim=dem_clim, alpha=0.6, **kwargs) + pltlib.add_cbar(axa[0,1], im, arr=src_dem_orig, clim=dem_clim, label=None) + axa[0,1].set_title('Source DEM') + + # just use stable surface here + arr = np.ones_like(stable_orig.data, dtype='uint8') + arr[stable_orig.mask == True] = 0 + # arr = np.ones_like(ref_dem_orig.data, dtype='uint8') + # arr[ref_dem_orig.mask == True] = 0 + # arr[src_dem_orig.mask == True] = 0 + axa[0,2].imshow(arr, clim=(0,1), cmap='gray', **kwargs) + axa[0,2].set_title('Surfaces for median shift') + + #This is empty + axa[0,3].axis('off') + + dz_clim = malib.calcperc_sym(diff_orig, (5, 95)) + im = axa[1,0].imshow(diff_orig, cmap='RdBu', clim=dz_clim, **kwargs) + pltlib.add_cbar(axa[1,0], im, arr=diff_orig, clim=dz_clim, label=None) + axa[1,0].set_title('Elev. Diff. Before (m)') + + im = axa[1,1].imshow(diff_align, cmap='RdBu', clim=dz_clim, **kwargs) + pltlib.add_cbar(axa[1,1], im, arr=diff_align, clim=dz_clim, label=None) + axa[1,1].set_title('Elev. Diff. After (m)') + + tight_dz_clim = (-2.0, 2.0) + im = axa[1,2].imshow(diff_align, cmap='RdBu', clim=tight_dz_clim, **kwargs) + pltlib.add_cbar(axa[1,2], im, arr=diff_align, clim=tight_dz_clim, label=None) + axa[1,2].set_title('Elev. Diff. After (m)') + + arr = np.ones_like(ref_dem_orig.data, dtype='uint8') + arr[ref_dem_orig.data == -9999] = 0 + updated_arr = np.ones_like(ref_dem_align.data, dtype='uint8') + updated_arr[ref_dem_align.data == -9999] = 0 + + bins = np.linspace(dz_clim[0], dz_clim[1], 128) + diff_orig_compressed = diff_orig[arr==1] + diff_align_compressed = diff_align[updated_arr==1] + axa[1,3].hist(diff_orig_compressed, bins, color='g', label='Before', alpha=0.5) + axa[1,3].hist(diff_align_compressed, bins, color='b', label='After', alpha=0.5) + axa[1,3].set_xlim(*dz_clim) + axa[1,3].axvline(0, color='k', linewidth=0.5, linestyle=':') + axa[1,3].set_xlabel('Elev. Diff. (m)') + axa[1,3].set_ylabel('Count (px)') + axa[1,3].set_title("Source - Reference") + + before_str = 'Before\nmed: %0.2f\nnmad: %0.2f' % (diff_orig_stats['med'], diff_orig_stats['nmad']) + axa[1,3].text(0.05, 0.95, before_str, va='top', color='g', transform=axa[1,3].transAxes, fontsize=8) + after_str = 'After\nmed: %0.2f\nnmad: %0.2f' % (diff_align_stats['med'], diff_align_stats['nmad']) + axa[1,3].text(0.65, 0.95, after_str, va='top', color='b', transform=axa[1,3].transAxes, fontsize=8); + + suptitle = '%s\nx: %+0.2fm, y: %+0.2fm, z: %+0.2fm' % (outprefix, dx, dy, dz) + f.suptitle(suptitle) + f.tight_layout() + plt.subplots_adjust(top=0.90) + + print(f"Writing out figure: {fig_fn}") + f.savefig(fig_fn, dpi=300) + +if __name__ == "__main__": + main() diff --git a/demcoreg/get_bareground.sh b/demcoreg/get_bareground.sh index fe97c70..22a564f 100755 --- a/demcoreg/get_bareground.sh +++ b/demcoreg/get_bareground.sh @@ -2,63 +2,113 @@ #David Shean #dshean@gmail.com +# Updated by JMH, 8.16.22 #This script will download global 30-m bare earth data #Needed for various masking applicaitons #Uses GNU parallel for processing -gdal_opt="-co TILED=YES -co COMPRESS=LZW -co BIGTIFF=YES" +# For more information, visit https://glad.umd.edu/dataset/global-2010-bare-ground-30-m +# Percentage of bare ground cover per output grid cell is encoded as integer values (1-100). +# Files are unsigned 8-bit data with spatial resolution of 0.00025° per pixel (approximately 30 meters per pixel at the equator). +# Data coverage is from 80° north to 60° south. It is divided into tiles of 10 x 10 degrees. +# The naming convention per tile refers to the latitude and longitude value of the upper left corner of the tile. +# Tiles over ocean area are included for completeness. + +# Updated links +# https://glad.umd.edu/Potapov/Bare_2010/[YY][N..S]_[XXX][E..W].tif +# Y ranges between 00–80 N +# Y ranges between 10–50 S +# X ranges between 010–180 W +# X ranges between 000–170 E +# increments by 10˚ in both lat and lon +gdal_opt="-co TILED=YES -co COMPRESS=LZW -co BIGTIFF=YES" if [ -z "$DATADIR" ] ; then - export DATADIR=$HOME/data + export DATADIR=$HOME/data/bare fi -cd $DATADIR +# Make datadir if it doesn't exist +if [ ! -e $DATADIR ] ; then + mkdir $DATADIR +fi -#~2010 global bare ground, 30 m -#When unzipped, this is 64 GB! -#http://landcover.usgs.gov/glc/BareGroundDescriptionAndDownloads.php +cd $DATADIR -be_zip_fn='bare2010.zip' -be_fn='bare2010' -be_vrt_fn='bare2010/bare2010.vrt' +# Final processed vrt of all tiles +be_vrt_fn=$DATADIR/bare2010.vrt -#This zipfile is 36 GB -if [ ! -e $be_fn ] ; then - url='http://edcintl.cr.usgs.gov/downloads/sciweb1/shared/gtc/downloads/bare2010.zip' - echo "Downloading $be_zip_fn" - wget -O $be_zip_fn $url +# If vrt DNE, download tiles with/without overwriting (default is no clobber) +overwrite=false +if [ ! -e $be_vrt_fn ] ; then + # This puts things into Potapov/Bare_2010 dir + # 51GB for 504 .tifs + echo "Downloading tiles" + if $overwrite ; then + time parallel -v wget -r -np -nH --no-check-certificate -A {} https://glad.umd.edu/Potapov/Bare_2010/ ::: "[0-9][0-9][NS]_[0-9][0-9][0-9][EW].tif" + else + time parallel -v wget -r -np -nH -nc --no-check-certificate -A {} https://glad.umd.edu/Potapov/Bare_2010/ ::: "[0-9][0-9][NS]_[0-9][0-9][0-9][EW].tif" + fi && + # move files to datadir + echo + echo "Moving tifs $DATADIR/Potapov/Bare_2010/*tif $DATADIR" + mv $DATADIR/Potapov/Bare_2010/*tif $DATADIR + echo + echo "Removing empty dirs" + rm -r $DATADIR/Potapov + #Now create a vrt for all input data + #Set nodata as 255, as 0 is valid bare ground percentage + echo "Building vrt of cleaned tiles: $be_vrt_fn" + gdalbuildvrt -vrtnodata 255 $be_vrt_fn $DATADIR/*.tif +else + echo + echo $be_vrt_fn exists fi -if [ ! -d $be_fn ] ; then - echo "Unzipping $be_zip_fn" - unzip $be_zip_fn -fi +# #~2010 global bare ground, 30 m +# #When unzipped, this is 64 GB! +# #http://landcover.usgs.gov/glc/BareGroundDescriptionAndDownloads.php -#These are 10x10-degree tiles at 30 m -#Original data are uncompressed, untiled -#Clean up -fn_list=$(ls bare2010/bare2010_v3/*bare2010_v3.tif) -c=$(wc -c $fn_list) -t=false -echo "Checking $c bare2010 files" -for fn in $fn_list -do - if ! gdalinfo $fn | grep -q LZW ; then - t=true - break - fi -done - -if $t ; then - echo "Cleaning up tiles" - parallel "gdal_translate $gdal_opt {} {.}_lzw.tif; mv {.}_lzw.tif {}" ::: bare2010/bare2010_v3/*bare2010_v3.tif -fi +# be_zip_fn='bare2010.zip' +# be_fn='bare2010' +# be_vrt_fn='bare2010/bare2010.vrt' -#Now create a vrt for all input data -#Set nodata as 255, as 0 is valid bare ground percentage -if [ ! -e $be_vrt_fn ] ; then - echo "Building vrt of cleaned tiles" - gdalbuildvrt -vrtnodata 255 $be_vrt_fn bare2010/bare2010_v3/*bare2010_v3.tif -fi +# #This zipfile is 36 GB +# if [ ! -e $be_fn ] ; then +# url='http://edcintl.cr.usgs.gov/downloads/sciweb1/shared/gtc/downloads/bare2010.zip' +# echo "Downloading $be_zip_fn" +# wget -O $be_zip_fn $url +# fi + +# if [ ! -d $be_fn ] ; then +# echo "Unzipping $be_zip_fn" +# unzip $be_zip_fn +# fi + +# #These are 10x10-degree tiles at 30 m +# #Original data are uncompressed, untiled +# #Clean up +# fn_list=$(ls $DATADIR/*.tif) +# c=$(wc -c $fn_list) +# t=false +# echo "Checking $c bare2010 files" +# for fn in $fn_list +# do +# if ! gdalinfo $fn | grep -q LZW ; then +# t=true +# break +# fi +# done + +# if $t ; then +# echo "Cleaning up tiles" +# parallel "gdal_translate $gdal_opt {} {.}_lzw.tif; mv {.}_lzw.tif {}" ::: $DATADIR/*.tif +# fi + +# #Now create a vrt for all input data +# #Set nodata as 255, as 0 is valid bare ground percentage +# if [ ! -e $be_vrt_fn ] ; then +# echo "Building vrt of cleaned tiles" +# gdalbuildvrt -vrtnodata 255 $be_vrt_fn $DATADIR/*.tif +# fi diff --git a/demcoreg/get_canopyheight.sh b/demcoreg/get_canopyheight.sh new file mode 100755 index 0000000..b1f3ab3 --- /dev/null +++ b/demcoreg/get_canopyheight.sh @@ -0,0 +1,101 @@ +#! /bin/bash + +# This script will download the global forest canopy height GEDI-Landsat data product detailed in https://doi.org/10.1016/j.rse.2020.112165 +# May be better than NLCD and global bareground data for identifying canopy cover and exposed rock surfaces + +# 7 regions available, see https://glad.umd.edu/dataset/gedi/ for more info and access to GEE app +# AUS: Australasia region [1.5 GB] +# NAFR: North Africa region [3.1 GB] (includes Europe and parts of Middle East) +# NAM: North America region [5.7 GB] +# NASIA: North Asia region [4.1 GB] +# SAFR: South Africa region [5.7 GB] +# SAM: South America region [6.9 GB] (includes part of Central America) +# SASIA: South Asia region [4.7 GB] + +# Usage +# ./get_canopyheight.sh [ region ] [ output data directory] + +################################################################################################### +# Help +Help() +{ + # Display Help + echo "Script downloader for 2019 Global Forest Canopy Height GEDI-Landsat 30 m products from" + echo "https://glad.umd.edu/dataset/gedi/" + echo + echo "Usage:" + echo "get_canopyheight.sh [ REGIONCODE ] [ output data directory ]" + echo + echo "Defaults" + echo " region NAM (North America)" + echo " datadir $HOME/data/canopyheight" + echo + echo "Available (continental) regions" + echo " AUS: Australasia region [1.5 GB]" + echo " NAFR: North Africa region [3.1 GB]" + echo " NAM: North America region [5.7 GB]" + echo " NASIA: North Asia region [4.1 GB]" + echo " SAFR: South Africa region [5.7 GB]" + echo " SAM: South America region [6.9 GB]" + echo " SASIA: South Asia region [4.7 GB]" + echo + echo "options:" + echo "h Print this Help." + echo +} +################################################################################################### + +# Get the options +while getopts ":h" option; do + case $option in + h) # display Help + Help + exit;; + esac +done + +################################################################################################### +set -e + +# Input region +region=$1 +DATADIR=$2 + +# Default is North America 'NAM' +if [ -z "$region" ] ; then + echo "Using default region: North America ['NAM']" + region='NAM' +fi + +gdal_opt="-co TILED=YES -co COMPRESS=LZW -co BIGTIFF=YES" + +if [ -z "$DATADIR" ] ; then + export DATADIR=$HOME/data/canopyheight +fi + +if [ ! -e $DATADIR ] ; then + mkdir $DATADIR +fi + +cd $DATADIR + +# Global 2019 grids, 30m from https://glad.geog.umd.edu/Potapov/Forest_height_2019/ +url="https://glad.geog.umd.edu/Potapov/Forest_height_2019/Forest_height_2019_${region}.tif" +canopyheight_fn="Forest_height_2019_${region}.tif" + +if [ ! -e $canopyheight_fn ] ; then + echo "Downloading $canopyheight_fn" + wget -O $canopyheight_fn $url +else + echo "Found existing ${canopyheight_fn}" +fi + +# Haven't tested this yet, checking to see if the bits are already as compressed as can be +if ! gdalinfo $canopyheight_fn | grep -q LZW ; then + echo "Compressing $canopyheight_fn" + gdal_translate -co TILED=YES -co COMPRESS=LZW -co BIGTIFF=IF_SAFER $canopyheight_fn temp.tif + if [ $? -eq 0 ] ; then + mv temp.tif $canopyheight_fn + fi +fi +################################################################################################### \ No newline at end of file diff --git a/demcoreg/get_nlcd.sh b/demcoreg/get_nlcd.sh index 480e5b3..b7ec883 100755 --- a/demcoreg/get_nlcd.sh +++ b/demcoreg/get_nlcd.sh @@ -5,17 +5,67 @@ #Should be better than global bareground data for identifying exposed rock surfaces # Usage: -# ./get_nlcd.sh year +# ./get_nlcd.sh [ region ] [ year ] [ ouput dat directory] + +################################################################################################### +# Help +Help() +{ + # Display Help + echo "Script downloader for NLCD 30 m products from" + echo "https://www.mrlc.gov/data" + echo + echo "Usage:" + echo "./get_nlcd.sh [ REGIONCODE ] [ year ] [ output data directory ]" + echo + echo "Defaults" + echo " region CONUS (Conterminous US)" + echo " yr 2019" + echo " datadir $HOME/data/" + echo + echo "Available regions" + echo " CONUS: Lower 48 US states" + echo " AK: Alaska" + echo + echo "options:" + echo "h Print this Help." + echo +} +################################################################################################### + +# Get the options +while getopts ":h" option; do + case $option in + h) # display Help + Help + exit;; + esac +done + +################################################################################################### set -e -# Input year -yr=$1 +echo ; echo Running: $0 $@ ; echo +# Input region +region=$1 +yr=$2 +DATADIR=$3 + +# Default is CONUS +if [ -z "$region" ] ; then + echo "Using default region: 'CONUS'" + region=CONUS +else + echo "Using input region: $region" +fi -# Default is 2016 +# Default is 2019 if [ -z "$yr" ] ; then - echo "Using default NLCD year 2016" - yr=2016 + echo "Using default NLCD year 2016" ; echo + yr=2019 +else + echo "Using input year: $yr" ; echo fi gdal_opt="-co TILED=YES -co COMPRESS=LZW -co BIGTIFF=YES" @@ -30,13 +80,21 @@ fi cd $DATADIR -#CONUS Land Cover (NLCD) grids, 30 m from https://www.mrlc.gov/data -url="https://s3-us-west-2.amazonaws.com/mrlc/nlcd_${yr}_land_cover_l48_20210604.zip" -nlcd_zip_fn="nlcd_${yr}_land_cover_l48_20210604.zip" -nlcd_fn="nlcd_${yr}_land_cover_l48_20210604.tif" +if [ $region == "CONUS" ] ; then + #CONUS Land Cover (NLCD) grids, 30 m from https://www.mrlc.gov/data + url="https://s3-us-west-2.amazonaws.com/mrlc/nlcd_${yr}_land_cover_l48_20210604.zip" +elif [ $region == "AK" ] ; then + #AK Land Cover (NLCD) grids, 30 m from https://www.mrlc.gov/data + url="https://s3-us-west-2.amazonaws.com/mrlc/NLCD_${yr}_Land_Cover_AK_20200724.zip" +else + echo "Region code not recognized: ${region}. Please enter 'CONUS' or 'AK'" +fi + +nlcd_zip_fn=$(basename ${url}) +nlcd_fn="${nlcd_zip_fn%.*}.tif" nlcd_dir=${nlcd_zip_fn%.*} -if [ ! -e $nlcd_zip_fn ] ; then +if [[ ! -f $nlcd_zip_fn && ! -f $nlcd_fn ]] ; then echo "Downloading $nlcd_zip_fn" wget -O $nlcd_zip_fn $url fi