Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
cmarshak opened this issue Jan 27, 2022 · 19 comments · Fixed by #47
Assignees
Labels
bug Something isn't working

Comments

@cmarshak
Copy link
Collaborator

cmarshak commented Jan 27, 2022

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:

  1. 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
  2. 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

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

More importantly, here is the data:

  1. 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
  2. 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

@cmarshak cmarshak added the bug Something isn't working label Jan 27, 2022
@sssangha
Copy link
Collaborator

sssangha commented Jan 27, 2022

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
Copy link
Collaborator Author

cmarshak commented Jan 27, 2022

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

@cmarshak
Copy link
Collaborator Author

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 cmarshak changed the title Incidence Angle Layer is not correctly read by ARIA-Tools Incidence Angle Layer is not correctly read by ARIA-Tools (specifically, Netcdf is not geo-referenced correctly) Jan 27, 2022
@cmarshak cmarshak changed the title Incidence Angle Layer is not correctly read by ARIA-Tools (specifically, Netcdf is not geo-referenced correctly) Incidence Angle Layer is not correctly read by ARIA-Tools [edit: Netcdf is not geo-referenced correctly] Jan 27, 2022
@cmarshak
Copy link
Collaborator Author

cmarshak commented Jan 28, 2022

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
Copy link
Collaborator Author

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

  1. Panoply correctly displays these datasets (clicking on the layer)
  2. We can view the netcdf layers (correctly georeferenced) in QGIS (as indicated in the screenshot below)
    image
  3. 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/*)

image

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

@cmarshak
Copy link
Collaborator Author

cmarshak commented Jan 31, 2022

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
Copy link
Collaborator Author

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
Copy link
Collaborator Author

cmarshak commented Jan 31, 2022

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
Copy link
Collaborator

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
Copy link
Collaborator Author

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
Copy link
Collaborator Author

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
Copy link
Collaborator Author

cmarshak commented Feb 1, 2022

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 cmarshak changed the title Incidence Angle Layer is not correctly read by ARIA-Tools [edit: Netcdf is not geo-referenced correctly] Incidence Angle Layer is not correctly read by ARIA-Tools [edit: buildVRT does not work with new GUNW products] Feb 1, 2022
@cmarshak
Copy link
Collaborator Author

cmarshak commented Feb 1, 2022

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
Copy link
Collaborator Author

cmarshak commented Feb 1, 2022

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
Copy link
Collaborator

dbekaert commented Feb 1, 2022 via email

@jhkennedy
Copy link
Collaborator

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
Copy link
Collaborator Author

cmarshak commented Feb 1, 2022

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

@cmarshak
Copy link
Collaborator Author

cmarshak commented Feb 1, 2022

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
Copy link
Collaborator Author

cmarshak commented Feb 2, 2022

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

5 participants