Incidence Angle Layer is not correctly read by ARIA-Tools [edit: buildVRT does not work with new GUNW products]

cmarshak commented 2 years ago

The issue currently is distilled in that gdal.buildVRT does not work on the /science/grids/imagingGeometry layers of the new GUNW product. Here is the minimum working example:

from pathlib import Path
test_gunw = Path('<gunw_id>.nc')
test_gunw_inc = f'NETCDF:"{test_gunw}":/science/grids/imagingGeometry/incidenceAngle'
test = gdal.BuildVRT('test.vrt', test_gunw_inc)
test = None

Trying to load test.vrt does not work (the xml file does correctly load the various layers). See this comment below for what the GUNW layers should look like.

Relevant data:

New data with incorrect georeferencing: https://hyp3-isce-contentbucket-4xpualmsjg98.s3.us-west-2.amazonaws.com/b4cb349d-b362-484f-ad06-2a35e1dec246/S1-GUNW-A-R-095-tops-20180705_20170920-045601-00170W_00052N-PP-ccb3-v2_0_5.nc
Old data with correct georeferencing (search result only): https://search.asf.alaska.edu/#/?searchType=List%20Search&searchList=S1-GUNW-A-R-137-tops-20210809_20210728-015824-35934N_33889N-PP-5ccf-v2_0_4&resultsLoaded=true&granule=S1-GUNW-A-R-137-tops-20210809_20210728-015824-35934N_33889N-PP-5ccf-v2_0_4-amplitude

Older introduction based on ARIA-tools (which needs `gdal.buildVRT`).

Based on the minimum working example which compares the previous GUNW products and the hyp3 version, this must be an issue with the current packing scripts in our Hyp3 Plugin (this repo).

This gist shows two minimum working examples:

https://gist.github.com/cmarshak/89c2826bf3334e89da6896d7ccf99f4e

It shows the new GUNWs do not allow "ariaExtract.py" for incidence angle
It shows the previous GUNWs allows for "ariaExtract.py" for incidence angle.

More importantly, here is the data:

New data with incorrect georeferencing: https://hyp3-isce-contentbucket-4xpualmsjg98.s3.us-west-2.amazonaws.com/b4cb349d-b362-484f-ad06-2a35e1dec246/S1-GUNW-A-R-095-tops-20180705_20170920-045601-00170W_00052N-PP-ccb3-v2_0_5.nc
Old data with correct georeferencing (search result only): https://search.asf.alaska.edu/#/?searchType=List%20Search&searchList=S1-GUNW-A-R-137-tops-20210809_20210728-015824-35934N_33889N-PP-5ccf-v2_0_4&resultsLoaded=true&granule=S1-GUNW-A-R-137-tops-20210809_20210728-015824-35934N_33889N-PP-5ccf-v2_0_4-amplitude

And the relevant command line script is:

<path_to_aria_tools>/tools/bin/ariaExtract.py --file 'S1-GUNW-A-R-095-tops-20180705_20170920-045601-00170W_00052N-PP-ccb3-v2_0_5.nc' --layer incidenceAngle -d <dem>.tif

Also, note that rasterio (same env) sees the older products as netcdf (correctly) but the newer versions as hdf5 (not correct); moreover, I can't load the transform of the new products. I was able to load the newer version in QGIS so maybe this is a rasterio issue. Still highlighting. This could partially be about rasterio, but maybe not.

The error is:

ARIA-tools Version: 1.1.0
*****************************************************************
*** Extract Product Function ***
*****************************************************************
All (1) GUNW products meet spatial bbox criteria.
Group GUNW products into spatiotemporally continuous interferograms.
All (1) interferograms are spatially continuous.
Thread count specified for gdal multiprocessing = 2
Applied cutline to produce 3 arc-sec SRTM DEM: ./DEM/SRTM_3arcsec.dem
Generating: incidenceAngle - [==================================================] 20180705_20170920Traceback (most recent call last):
  File "/home/cmarshak/leffe2-cmarshak/ARIA-tools/tools/bin/ariaExtract.py", line 21, in <module>
    main(inps)
  File "/mnt/leffe-data2/cmarshak/ARIA-tools/tools/ARIAtools/extractProduct.py", line 1196, in main
    export_products(standardproduct_info.products[1], standardproduct_info.bbox_file,
  File "/mnt/leffe-data2/cmarshak/ARIA-tools/tools/ARIAtools/extractProduct.py", line 622, in export_products
    finalize_metadata(outname, bounds, dem_bounds, prods_TOTbbox, dem, lat, lon, mask, outputFormat, verbose=verbose)
  File "/mnt/leffe-data2/cmarshak/ARIA-tools/tools/ARIAtools/extractProduct.py", line 698, in finalize_metadata
    data_array = metadata_qualitycheck(data_array, os.path.basename(os.path.dirname(outname)), outname, verbose).data_array
  File "/mnt/leffe-data2/cmarshak/ARIA-tools/tools/ARIAtools/extractProduct.py", line 132, in __init__
    self.__run__()
  File "/mnt/leffe-data2/cmarshak/ARIA-tools/tools/ARIAtools/extractProduct.py", line 233, in __run__
    if min(rsquaredarr_rng) < 0.97 and max(std_errarr_rng) > 0.0015:
ValueError: min() arg is an empty sequence

or here in the ARIA-Tools: https://github.com/aria-tools/ARIA-tools/blob/dev/tools/ARIAtools/extractProduct.py#L149-L220

sssangha commented 2 years ago

This leads to an issue in ARIA-tools as the initial gdal.BuildVRT call that creates a VRT pointing to the metadata layers doesn't work as expected.

E.g. with the newer products if you run the following

from osgeo import gdal

fname = 'NETCDF:"/u/leffe-data/ssangha/Charlie_test/test2/products/S1-GUNW-A-R-095-tops-20191016_20190910-045612-00170W_00052N-PP-fbd7-v2_0_5.nc":/science/grids/imagingGeometry/incidenceAngle' gdal.BuildVRT('test.vrt', fname)

The VRT that is generated doesn't contain relative paths:

<VRTDataset rasterXSize="49" rasterYSize="28">
  <SRS dataAxisToSRSAxisMapping="2,1">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG"
,"9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS>
  <GeoTransform> -1.6504999999999998e+02, -1.0000000000000024e-01,  0.0000000000000000e+00,  5.4950000000000003e+01,  0.0000000000000000e+00, -1.0000000000000010e-01</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>0</NoDataValue>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="2">
    <NoDataValue>0</NoDataValue>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="3">
    <NoDataValue>0</NoDataValue>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="4">
    <NoDataValue>0</NoDataValue>
  </VRTRasterBand>
</VRTDataset>

As opposed to the older products which generate the expected VRT and don't lead to issues in metadata layer extraction:

<VRTDataset rasterXSize="37" rasterYSize="27">
 <SRS dataAxisToSRSAxisMapping="2,1">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG"
,"9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS>
 <GeoTransform> -1.2145000000000000e+02, 9.9999999999999839e-02, 0.0000000000000000e+00, 3.6250000000000000e+01, 0.0000000000000000e+00, -1.0000000000000006e-01</GeoTransform>
 <VRTRasterBand dataType="Float32" band="1">
  <NoDataValue>0</NoDataValue>
  <ComplexSource>
   <SourceFilename relativeToVRT="0">NETCDF:"/u/leffe-data/ssangha/Charlie_test/test2/CA_test/products/S1-GUNW-A-R-137-tops-20150811_20150507-015829-35934N_33888N-PP-472f-v2_0_2.nc":/science/grids/imagingGeometry/incidenceAngle</SourceFilename>
   <SourceBand>1</SourceBand>
   <SourceProperties RasterXSize="37" RasterYSize="27" DataType="Float32" BlockXSize="37" BlockYSize="27" />
   <SrcRect xOff="0" yOff="0" xSize="37" ySize="27" />
   <DstRect xOff="0" yOff="0" xSize="37" ySize="27" />
   <NODATA>0</NODATA>
  </ComplexSource>
 </VRTRasterBand>
 <VRTRasterBand dataType="Float32" band="2">
  <NoDataValue>0</NoDataValue>
  <ComplexSource>
   <SourceFilename relativeToVRT="0">NETCDF:"/u/leffe-data/ssangha/Charlie_test/test2/CA_test/products/S1-GUNW-A-R-137-tops-20150811_20150507-015829-35934N_33888N-PP-472f-v2_0_2.nc":/science/grids/imagingGeometry/incidenceAngle</SourceFilename>
   <SourceBand>2</SourceBand>
   <SourceProperties RasterXSize="37" RasterYSize="27" DataType="Float32" BlockXSize="37" BlockYSize="27" />
   <SrcRect xOff="0" yOff="0" xSize="37" ySize="27" />
   <DstRect xOff="0" yOff="0" xSize="37" ySize="27" />
   <NODATA>0</NODATA>
  </ComplexSource>
 </VRTRasterBand>
 <VRTRasterBand dataType="Float32" band="3">
  <NoDataValue>0</NoDataValue>
  <ComplexSource>
   <SourceFilename relativeToVRT="0">NETCDF:"/u/leffe-data/ssangha/Charlie_test/test2/CA_test/products/S1-GUNW-A-R-137-tops-20150811_20150507-015829-35934N_33888N-PP-472f-v2_0_2.nc":/science/grids/imagingGeometry/incidenceAngle</SourceFilename>
   <SourceBand>3</SourceBand>
   <SourceProperties RasterXSize="37" RasterYSize="27" DataType="Float32" BlockXSize="37" BlockYSize="27" />
   <SrcRect xOff="0" yOff="0" xSize="37" ySize="27" />
   <DstRect xOff="0" yOff="0" xSize="37" ySize="27" />
   <NODATA>0</NODATA>
  </ComplexSource>
 </VRTRasterBand>
 <VRTRasterBand dataType="Float32" band="4">
  <NoDataValue>0</NoDataValue>
  <ComplexSource>
   <SourceFilename relativeToVRT="0">NETCDF:"/u/leffe-data/ssangha/Charlie_test/test2/CA_test/products/S1-GUNW-A-R-137-tops-20150811_20150507-015829-35934N_33888N-PP-472f-v2_0_2.nc":/science/grids/imagingGeometry/incidenceAngle</SourceFilename>
   <SourceBand>4</SourceBand>
   <SourceProperties RasterXSize="37" RasterYSize="27" DataType="Float32" BlockXSize="37" BlockYSize="27" />
   <SrcRect xOff="0" yOff="0" xSize="37" ySize="27" />
   <DstRect xOff="0" yOff="0" xSize="37" ySize="27" />
   <NODATA>0</NODATA>
  </ComplexSource>
 </VRTRasterBand>
</VRTDataset>

cmarshak commented 2 years ago

I think its related to #44 . Darn...

cmarshak commented 2 years ago

And finally, here is the dockerfile for the older topsApp "plugin" - which installs netcdf library via pip.

https://github.com/aria-jpl/topsApp_pge/blob/develop/docker/Dockerfile

cmarshak commented 2 years ago

Here are some logs using this gist to demonstrate issues with packaging.

make_geocube.log packaging.log

There are pretty standard pyproj warnings.

This illustrates the packaging workflow.

There are some additional raster layers generated by the CL script makeGeocube.py and then they are packaged together using a different command line script. The former rasters are included in a metadata.h5 file in the processing directory.

cmarshak commented 2 years ago

To summarize, it's hard to say if this netcdf issue persisted throughout the development of this plugin because:

Panoply correctly displays these datasets (clicking on the layer)
We can view the netcdf layers (correctly georeferenced) in QGIS (as indicated in the screenshot below)
ariaExtract (not specifying the incidence angle) works as we expect.

I am not sure of the best way to proceed other than to preserve the existing packaging and manually update the relevant raster layers using xarray to ensure correct georeferencing. The relevant layers in the screenshot (i.e./science/grids/data/* and /science/grids/imagingGeometry/*)

If there is some other environment issue that could be at play, it would great to now a roadmap.

cmarshak commented 2 years ago

So, I installed rioxarray and xarray and was able to correctly read the geotransform with rasterio - so I am a little lost as to what is wrong in the product again; what precisely is missing for this to work as we expect it to?

In [1]: import rasterio

with rasterio.open('S1-GUNW-A-R-095-tops-20200823_20191004-045617-00170W_00052N-PP-6e20-v2_0_5.nc') as ds:
    subdatasets = ds.subdatasets

with rasterio.open(subdatasets[2]) as ds:
    t = ds.transform

t * (0, 0)
Out[1]: (-169.678625, 54.6120833305)

To be clear, this was a change on probing the product, not generating it. Still, the error with ARIA-tools reported in the beginning of this issue.

cmarshak commented 2 years ago

The buildvrt issue remains and maybe it is isolated there. When I call info on the newer and older products they seem identical.

https://gist.github.com/cmarshak/5948c84962c50c4a37e8e8fce358214d

cmarshak commented 2 years ago

My guess is that it has something to do with pyproj because the latitude and longitudes are switched in the Info command. See the gist in the previous comment.

Newer version:

'coordinateSystem': {'wkt': 'GEOGCRS["WGS 84",\n    DATUM["World Geodetic System 1984",\n        ELLIPSOID["WGS 84",6378137,298.257223563,\n            LENGTHUNIT["metre",1]]],\n    PRIMEM["Greenwich",0,\n        ANGLEUNIT["degree",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS["latitude",north,\n            ORDER[1],\n            ANGLEUNIT["degree",0.0174532925199433]],\n        AXIS["longitude",east,\n            ORDER[2],\n            ANGLEUNIT["degree",0.0174532925199433]],\n    ID["EPSG",4326]]' ...

vs.

Older version:

'coordinateSystem': {'wkt': 'GEOGCRS["WGS 84",\n    DATUM["World Geodetic System 1984",\n        ELLIPSOID["WGS 84",6378137,298.257223563,\n            LENGTHUNIT["metre",1]]],\n    PRIMEM["Greenwich",0,\n        ANGLEUNIT["degree",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS["longitude",east,\n            ORDER[1],\n            ANGLEUNIT["degree",0.0174532925199433]],\n        AXIS["latitude",north,\n            O ...

jhkennedy commented 2 years ago

This is due to an upgrade of the pyproj library, and a "fun" gotacha:

https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6

Fix is to go through the packaging_utils and update any proj/transformation calls to the new interface, and/or the same in AIRA-tools

cmarshak commented 2 years ago

Ok, David brought this to my attention when we updated the packaing scripts for the newest ISCE2 (e.g. 'master' --> 'reference'). Looks like this issue was "sneaking" by to the final products because the product worked in the ways we had outlined to inspect the product (e.g. QGIS and panoply as seen above).

I am confident it is occurring in the makeGeocube.py https://github.com/ACCESS-Cloud-Based-InSAR/DockerizedTopsApp/blob/3e9aea8580f98b33c110410c418f9dc8e5124c7a/isce2_topsapp/packaging_utils/makeGeocube.py

I will update this and test out the new product.

cmarshak commented 2 years ago

The VRT issue that @sssangha still does not work even after updating the pyproj and removing those warnings from makeGeocube.py. I am not sure what the issue is here.

cmarshak commented 2 years ago

To determine gdal issue that is causing gdal.buildVRT to fail, we are submitting a new GUNW job for the older version highlighted above to compare layer by layer.

Specifically,

ref: S1B_IW_SLC__1SDV_20210809T015810_20210809T015838_028163_035C1C_1C8D

sec: S1B_IW_SLC__1SDV_20210728T015742_20210728T015812_027988_0356CA_5433 S1B_IW_SLC__1SDV_20210728T015835_20210728T015902_027988_0356CA_D3EB

again for this GUNW (2.0.4): https://search.asf.alaska.edu/#/?searchType=List%20Search&searchList=S1-GUNW-A-R-137-top[…]09_20210728-015824-35934N_33889N-PP-5ccf-v2_0_4-amplitude

cmarshak commented 2 years ago

I believe that I found the source of this issue.

The ImagingGeometry rasters have a non-standard geotransform.

import rasterio

with rasterio.open('S1-GUNW-A-R-095-tops-20200823_20191004-045617-00170W_00052N-PP-6e20-v2_0_5.nc') as ds:
    subdatasets = ds.subdatasets

with rasterio.open(subdatasets[-2]) as ds:
    t = ds.transform
    p = ds.profile
    X = ds.read()

t.to_gdal()
# (-165.04999999999998,  -0.10000000000000024, 0.0, 54.95,  0.0, -0.1000000000000001)

That second number should be positive for most gdal geotransforms. Here is a standard one from the other product (-121.45, 0.09999999999999984, 0.0, 36.25, 0.0, -0.10000000000000006).

I hope we can figure out where this is getting written.

Relevant pieces of the plugin can roughly be determined from this template:

https://github.com/ACCESS-Cloud-Based-InSAR/DockerizedTopsApp/blob/dev/isce2_topsapp/templates/nc_packaging_template.json#L141-L251

These functions are then in https://github.com/ACCESS-Cloud-Based-InSAR/DockerizedTopsApp/blob/dev/isce2_topsapp/packaging_utils/isce_functions.py

cmarshak commented 2 years ago

I have a draft fix for this bug. It's ugly and honestly don't know how it was introduced (by me).

Hoping its correct now.

We will clearly have to do more tests to assure its correct. May want to make it slightly prettier.

See this: https://github.com/ACCESS-Cloud-Based-InSAR/DockerizedTopsApp/pull/47/files#diff-f99602fd7669b76bfe8bed3563ddcbb48a34ae8c3acb404584b93202f9c7ea27

dbekaert commented 2 years ago

Yes believe you are correct. If i recall geotrans in gdal when requested should be consistent top left corner of a product. So second pos and 5th neg would be north south east west stored image.

@ sim can you also confirm grid spacings to be consistent. Should be integer increment with dem spacing of old products

On Mon, Jan 31, 2022 at 6:47 PM Charlie Marshak @.***> wrote:

I have a draft fix for this bug. It's ugly and honestly don't know how it was introduced (by me).

Hoping its correct now.

We will clearly have to do more tests to assure its correct. May want to make it slightly prettier.

— Reply to this email directly, view it on GitHub https://github.com/ACCESS-Cloud-Based-InSAR/DockerizedTopsApp/issues/46#issuecomment-1026426383, or unsubscribe https://github.com/notifications/unsubscribe-auth/AESZPSKMPATZLBLQVFXTLY3UY5CUZANCNFSM5M7DADHA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you are subscribed to this thread.Message ID: <ACCESS-Cloud-Based-InSAR/DockerizedTopsApp/issues/46/1026426383@ github.com>

jhkennedy commented 2 years ago

The example secondary scenes are missing a granule resulting in a discontinuous secondary frame (also missing from v204 GUNW metadata). Should be:

reference scenes:

S1B_IW_SLC__1SDV_20210809T015810_20210809T015838_028163_035C1C_1C8D

secondary scenes:

S1B_IW_SLC__1SDV_20210728T015742_20210728T015812_027988_0356CA_5433
S1B_IW_SLC__1SDV_20210728T015809_20210728T015837_027988_0356CA_917C
S1B_IW_SLC__1SDV_20210728T015835_20210728T015902_027988_0356CA_D3EB

cmarshak commented 2 years ago

We will make sure with some sample products this works correctly.

cmarshak commented 2 years ago

I re-ran the newest dev branch from end-to-end using this notebook: https://gist.github.com/cmarshak/77c251a02c8985153a1c8110b9d25bfe

Confident the issue has been resolved.

For further inspection, here is the same GUNW 2.0.5 from hyp3 deployment: https://hyp3-isce.asf.alaska.edu/jobs/f07519ee-c239-4b4c-97d0-cd49b0f60566

cmarshak commented 2 years ago

From the above single job - I successfully re-ran the above analysis and verified both ARIA-tools works and that the comparisons with the previous version also works.

We can re-open down the road if errors were not properly addressed.

ACCESS-Cloud-Based-InSAR / DockerizedTopsApp