Large file sizes and missing 'Elevation' layer name

[jlwoodr]

Hello,

I am trying to convert a large set of elevation raster tiles that I bulk downloaded from NOAA Digital Coast. Unfortunately, these tiles do not label their band names so vyperdatum fails to convert them from NAVD88 to MSL. I created a janky work around by just adding '' to the list of acceptable raster elevation layer names, but unfortunately the output is a new raster with 2 bands and over double the file size. I simply want to take my 1 banded elevation raster and convert it to LMSL. Is there a way to do this? Thank you for your help!

Johnathan

[ericgyounkin]

Hi @jlwoodr,

I talked to my coworker who is more involved with the raster side of things. It sounds like we just need to add an optional argument to allow you to specify the band name. You'd still need to actually know which band name corresponds with depth/elevation for the raster. Would that work?

[jlwoodr]

Hey Eric,

Thank you for your response. When I look at the bands in QGIS the below screen shot is what is looks like. It just says 'Band 1', not 'Band 1 : Elevation'. I can also send an email with the raster file I am looking at in case ya'll would like to take a look. I do not think I can attach it here though unfortunately. Thank you for your help!

[jlwoodr]

For reference, the vyperdatum test raster shows the following:

[ericgyounkin]

Can you go to Raster - Miscellaneous - Raster Information, Run the GDAL info command, and tell me what you see for your Band 1?

Should look something like this (you'll have to scroll down)

Band 1 Block=128x8 Type=Float32, ColorInterp=Gray Description = Depth Min=52.005 Max=70.657 Minimum=52.005, Maximum=70.657, Mean=56.015, StdDev=4.071 NoData Value=1000000 Unit Type: m Metadata: STATISTICS_MAXIMUM=70.656997680664 STATISTICS_MEAN=56.015299130917 STATISTICS_MINIMUM=52.005001068115 STATISTICS_STDDEV=4.0713242080392

We are looking for the Description tag.

[jlwoodr]

Here is the printout from that command

Driver: GTiff/GeoTIFF Files: C:/Users/jlwoodr3/Desktop/Oceanmesh_data_and_scripts/NC/vdatum_test/ncei19_n34x75_w076x50_2019v2.tif Size is 8112, 8112 Coordinate System is: GEOGCRS["NAD83", DATUM["North American Datum 1983", ELLIPSOID["GRS 1980",6378137,298.257222101004, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4269]] Data axis to CRS axis mapping: 2,1 Origin = (-76.500185185199996,34.750185185200003) Pixel Size = (0.000030864197535,-0.000030864197535) Metadata: AREA_OR_POINT=Area TIFFTAG_COPYRIGHT=DOC/NOAA/NESDIS/NCEI > National Centers for Environmental Information, NESDIS, NOAA, U.S. Department of Commerce TIFFTAG_DATETIME=2019:04:23 00:00:00 TIFFTAG_IMAGEDESCRIPTION=Topography-Bathymetry; NAVD88 Image Structure Metadata: COMPRESSION=DEFLATE INTERLEAVE=BAND PREDICTOR=3 Corner Coordinates: Upper Left ( -76.5001852, 34.7501852) ( 76d30' 0.67"W, 34d45' 0.67"N) Lower Left ( -76.5001852, 34.4998148) ( 76d30' 0.67"W, 34d29'59.33"N) Upper Right ( -76.2498148, 34.7501852) ( 76d14'59.33"W, 34d45' 0.67"N) Lower Right ( -76.2498148, 34.4998148) ( 76d14'59.33"W, 34d29'59.33"N) Center ( -76.3750000, 34.6250000) ( 76d22'30.00"W, 34d37'30.00"N) Band 1 Block=512x512 Type=Float32, ColorInterp=Gray Min=-36.911 Max=5.204 Minimum=-36.911, Maximum=5.204, Mean=-19.715, StdDev=7.912 NoData Value=-9999 Overviews: 4056x4056, 2028x2028, 1014x1014, 507x507, 254x254 Unit Type: metre Metadata: STATISTICS_MAXIMUM=5.204381942749 STATISTICS_MEAN=-19.715157438006 STATISTICS_MINIMUM=-36.911239624023 STATISTICS_STDDEV=7.9124649044338

[jlwoodr]

It appears that there is no "description" tag for these NCEI rasters.

[ericgyounkin]

Yes, like the band name is just "" like you discovered. I think that probably still works for my idea of just giving an override_band_name option, let me look into it a bit more.

The other thing I'm seeing from your initial post is that you only want the resultant depth, not the uncertainty/contributor layers. That is not currently an option. It will return all three layers regardless.

[jlwoodr]

Eric, I appreciate your help. Yeah I am trying to translate over 400 data sets so I am really trying to limit the size of my files. Maybe I can add a post processing function to rewrite the raster with only one band?

[ericgyounkin]

I might be able to add something there to support. Let me look into it right after I get this other thing sorted out.

[jlwoodr]

Eric, Sounds great! No rush!

[ericgyounkin]

I will also say that if you are just wanting to just process files and are not interested in the uncertainty/CRS stuff that we built into vyperdatum, you could also try the VDatum application. They added raster file processing fairly recently. I haven't tested it, but it might be a quick fix for you.

You just run the vdatum.jar application in the downloaded vdatum folder.

[ericgyounkin]

Hello,

I just pushed an update for 0.1.12, you can get it by pip installing the repo again.

https://github.com/noaa-ocs-hydrography/vyperdatum/commit/f2e2d2c56ba828cb7a878d38a68a3d7a4978dc2e

I've added logic for assuming that the band in a one band raster is depth/elevation. I think that is going to always be true.

I've also added the ability to skip uncertainty in transform. I ran it on an example dataset that I have, and it appeared to work fine. See below.

from vyperdatum.raster import VyperRaster

vr = VyperRaster(r"C:\collab\dasktest\data_dir\outputtest\tj_patch_test_710_20220416_134440\tj_patch_test_710_20220104_203429.tif")
Operating on C:\collab\dasktest\data_dir\outputtest\tj_patch_test_710_20220416_134440\tj_patch_test_710_20220104_203429.tif

vr.transform_raster((6347, 'tss'), (6347, 'navd88'), output_filename=r'C:\collab\dasktest\data_dir\outputtest\tj_patch_test_710_20220416_134440\tj_patch_test_710_20220104_203429_tss.tif', include_uncertainty=False)
Begin work on tj_patch_test_710_20220104_203429.tif
Unable to find uncertainty or vertical uncertainty layer by name, layers=['depth']
Unable to find contributor layer by name, layers=['depth']
Applying vdatum separation model to 32768 total points
Raster transformation complete: Elapsed time 0.23219359999999867 seconds
Out[4]: 
((array([[1000000., 1000000., 1000000., ..., 1000000., 1000000., 1000000.],
         [1000000., 1000000., 1000000., ..., 1000000., 1000000., 1000000.],
         [1000000., 1000000., 1000000., ..., 1000000., 1000000., 1000000.],
         ...,
         [1000000., 1000000., 1000000., ..., 1000000., 1000000., 1000000.],
         [1000000., 1000000., 1000000., ..., 1000000., 1000000., 1000000.],
         [1000000., 1000000., 1000000., ..., 1000000., 1000000., 1000000.]]),
  None,
  None),
 ['Depth'],
 [1000000.0])

[jlwoodr]

Eric,

You're the man, thank you! I just get this one error at the end.

[ericgyounkin]

Ah I see. Fixed this in 0.1.13. I missed an 'include_uncertainty' stanza in the allow points outside coverage routine.

Which I guess means you have some points outside VDatum coverage, not sure if that is expected for your dataset?

[jlwoodr]

Eric,

Okay that seemed to work. Thanks! I also have another question about odd numbers in the output dataset. For example, in the test dataset provided on github, I am getting areas in band 1 with radically different elevation (in the input an area might have -17 m elevation and in the output it is 8 m). What could be causing this? I think I am just unsure of what vdatum is exactly doing.

[ericgyounkin]

Vyperdatum is going to shift the input raster vertically from the datum you provide to the datum you want. What are you providing for input/output datum in transform_raster?

[jlwoodr]

For input datum I am providing NAVD88 and output LMSL ('tss') in vyperdatum I believe.

[ericgyounkin]

Hi @jlwoodr,

First, recommend you update again, I had to fix one last thing with the layer indexing. 0.1.14.

I'm a little unused the raster side of vyperdatum. Let me put up a comprehensive example here to illustrate.

Can you try running this as shown below? I found that I got a strange result when I ran it like I showed you before, including the EPSG code in the datum definition. Providing just 'tss' and 'navd88' seemed to make it work ok.

from vyperdatum.raster import VyperRaster

vr = VyperRaster(r"C:\Pydro22_Dev\NOAA\site-packages\Python38\git_repos\vyperdatum\data\tiff\test.tiff")
Operating on C:\Pydro22_Dev\NOAA\site-packages\Python38\git_repos\vyperdatum\data\tiff\test.tiff

vr.transform_raster('tss', 'navd88', output_filename=r"C:\Pydro22_Dev\NOAA\site-packages\Python38\git_repos\vyperdatum\data\tiff\test_tss.tiff", include_uncertainty=False)

End result is a shift of about 10 cm or so. Which is basically the vdatum answer for tss-navd88.

[jlwoodr]

Eric,

This fix seems to work great! I really appreciate all of your help, thanks a ton!