Canopy tiff is relative to snow depth instead of snow -off

[ZachHoppinen]

@tatemeehan - reports that the canopy tiff is relative to snow raster rather than snow-off reference dem. @brentwilder or @shadoneel is this something either of you have noticed?

In https://github.com/SnowEx/ice-road-copters/blob/main/scripts/ice-road-pipeline.py lines 129-135

# difference two rasters to find snow depth
ref_dem_path = join(results_dir, 'dem.tif')
canopy_fp = join(ice_dir, f'{basename(in_dir)}-canopyheight.tif')
canopy = rio.open_rasterio(canopy_tif, masked=True) 
matched = canopy.rio.reproject_match(snowoff)
canopyheight = matched - snowoff
canopyheight.rio.to_raster(canopy_fp)

Seems like the pipeline should be right so I am wondering if the snowoff variable is getting overwritten somewhere in the scripts?

[ZachHoppinen]

lines 121 and 123:

 ref_dem_path = join(results_dir, 'dem.tif')
snowoff = rio.open_rasterio(ref_dem_path, masked=True)

So maybe when we copy in the reference dem to dem.tif in the results directory it is getting changed to snow depth. We should be able to check one of the results directories to see if the dem.tif is snow raster or snow-off raster.

[ZachHoppinen]

March 17th, 2022 for mores is the site he saw it in.

[ZachHoppinen]

Just pulled the file from: kelvin:/SNOWDATA/IDALS/REF_DEM/MCS_REFDEM_WGS84.tif and kelvin:/SNOWDATA/IDALS/2022/20220317_MCS/ice-road/results/dem.tif and kelvin:/SNOWDATA/IDALS/2022/20220317_MCS/ice-road/20220317_MCS-snow.tif

and they are all identical DEMs. Weird

[ZachHoppinen]

So the problem is that the first returns are from the snow surface once the vegetation is below the snow surface. Two options:

1. Mask points where veg height == snow height to avoid giving the impression we know the snow height
2. subtract the snow depth from the veg height to make those areas 0.

[brentwilder]

@ZachKeskinen @tatemeehan @shadoneel

Thanks for catching this Tate. We kind of through in the canopy product at the end as a bonus to the snow depth product, with the intention that canopy would need further work/testing. This gdal_calc line shown here should do what we need.

# difference two rasters to find **canopy height**
ref_dem_path = join(results_dir, 'dem.tif')
canopy_fp = join(ice_dir, f'{basename(in_dir)}-canopyheight.tif')
canopy = rio.open_rasterio(canopy_tif, masked=True) 
matched = canopy.rio.reproject_match(snowoff)
canopyheight = matched - snowoff
canopyheight.rio.to_raster(canopy_fp)

# Apply a threshold.. was finding that setting it to zero was not quite working because small amounts of noise
# But I was able to get a clean CHM by adding 0.1m to the threshold.
canopy_clean_fp = join(ice_dir, f'{basename(in_dir)}-canopyheight_clean.tif')
os.system(f'gdal_calc.py -A {canopy_fp} -B {snow_depth_path} \
                --outfile={canopy_clean_fp} --type Float64 --overwrite \
                --calc="where((A<=B+0.1),0,A)" --quiet')

[brentwilder]

Before threshold:

After threshold:

[ZachHoppinen]

Nice Brent! That looks great, and I think we should roll that into the code if you want to make a branch and pull it in.

We should go with yours since you have tested it but the xarray version would look like:

canopyheight = matched - snowoff

# mask snow depth from vegetation
canopyheight = canopyheight.where(canopyheight < snowon + 0.1)
# might not need this but good to make sure masking didn't add anything weird to the nans
canopyheight = canopyheight.where(~snowon.isnull())

canopyheight.rio.to_raster(canopy_fp)

[brentwilder]

@ZachKeskinen I like trying to keep it in rasterio. I believe the sign is flipped - should read canopyheight > snowon + 0.1. Slight problem is that this then creates NaN in the dataset where we should have zeros for the CHM. These zeros are burned in with the gdal logical operation, but I don't know how to do this in rasterio? Do you have any ideas?

**EDIT: Actually , I like this better.. Looking at this CHM. It's not super complete since this wasn't the intention. For example, there are data missing in the center of large tree canopies. I think leaving as NaN would be sufficient and perhaps better even.

[brentwilder]

@ZachKeskinen I don't have permissions to Kelvin code but I will push it locally

[ZachHoppinen]

Awesome, great catch on that sign Brent. I agree that NANs are probably a better option than 0s since there could be vegetation that we just don't know the height of.

I think Shad is who cloned that repo on /SNOWDATA/IDALS so @shadoneel you will need to pull in any code changes.

[brentwilder]

Wait.. the center of the canopies are becoming NaN because the snow depth there is NaN.. and the "where" logical statement can't compute those areas. we need another line to catch these. @ZachKeskinen

[brentwilder]

Got it!

# mask snow depth from vegetation
canopyheight = canopyheight.where((canopyheight > snowon + 0.1) | (snowon.isnull()))

[brentwilder]

[ZachHoppinen]

Beautiful!

[shadoneel]

Hey guys I’m in a long set of meetings this week but hate to slow progress. Can you tell me exactly what needs doing and I can try tomorrow? Sent from my iPhoneOn Mar 6, 2023, at 11:48 AM, Zachary Keskinen @.***> wrote: Awesome, great catch on that sign Brent. I agree that NANs are probably a better option than 0s since there could be vegetation that we just don't know the height of. I think Shad is who cloned that repo on /SNOWDATA/IDALS so @shadoneel you will need to pull in any code changes.

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: @.***>

[brentwilder]

@shadoneel if using VS Code, you just need to click "Sync Changes"