Help needed to Geocode the dolphin output

[Shima-shoki]

The answer to this question might have already been answered by #544 , and I referenced to the material, but I'm not sure how to apply to dolphin. Please kindly help me if possible to resolve this step, as I want to process non-OPERA CSLC dataset outside of the north America, so I need to do geocoding by myself.

I have run dolphin to get the final velocity map. My guess is that we would be able to use the lat.rdr.full.vrt/lon.rdr.full.vrt inside the geom_reference folder, which was produced by stackSentinel.py. However, it seems that the shape of the velocity map and the lon.rdr.full.vrt are not the same. On the other hand, the shape of the co-registered data is consistent with the lon.rdr.full.vrt. Would it be possible to still use the longitude (latitude) map for georeferencing the final dolphin output, or should I (somehow) geocode the SLC data before the dolphin application? Any help is very much appreciated 🙏

[scottstanie]

The answer to this question might have already been answered by #544 , and I referenced to the material, but I'm not sure how to apply to dolphin. Please kindly help me if possible to resolve this step, as I want to process non-OPERA CSLC dataset outside of the north America, so I need to do geocoding by myself.

I have run dolphin to get the final velocity map. My guess is that we would be able to use the lat.rdr.full.vrt/lon.rdr.full.vrt inside the geom_reference folder, which was produced by stackSentinel.py. However, it seems that the shape of the velocity map and the lon.rdr.full.vrt are not the same. On the other hand, the shape of the co-registered data is consistent with the lon.rdr.full.vrt. Would it be possible to still use the longitude (latitude) map for georeferencing the final dolphin output, or should I (somehow) geocode the SLC data before the dolphin application? Any help is very much appreciated 🙏

It should definitely be okay to run dolphin in radar coordinates.

What was the configuration you used for this? Did you set "strides" to be larger than 1?

[scottstanie]

You could probably still geocode the results with the strided data. I think you can just do something like

gdal_translate -tr 3 3 lon.rdr.full.vrt lon.strided.vrt
gdal_translate -tr 3 3 lat.rdr.full.vrt lat.strided.vrt

I dont have one of these ISCE2 results handy to ensure that GDAL reads lon.rdr.full.vrt as having 1.0 resolution. But the idea is to make a downsampled version of lat and lon which matches your velocity, and then run a export / geocode script like the UNAVCO notebook.

[Shima-shoki]

I see, thank you very much for the interesting solution!
The stride 1 looks taking a huge GPU memory space, so if the larger strides work, it's very convenient 👍
I'll first try with the stride 1 & lon/lat rdr.full.vrt to geocode, and if they're successful, next step would be the geocoding with the stride options :)

[scottstanie]

Yes I mostly process with strides > 1 for that reason. But on that note (and something that could go on the future FAQ) - if you do strides of 1, you can still make the RAM usage lower by lowering the block_shape (within worker_settings) to something like (256, 256), or even lower if you have a very small GPU.

[Shima-shoki]

@scottstanie
Based on the comments, I was able to geocode the dolphin output with strides👍
The "-tr" option in gdal_translate didn't work, unfortunately, because the dataset doesn't have proper geometric properties. However, the output file shape can be specified with -outsize option, and the things went well afterwards!

#Dolphin input command

!dolphin config --slc-files /content/drive/MyDrive/SatelliteImageProcessings/InSAR_Python/ISCE2-GPU_SLCstack/merged/SLC/*/*.vrt --strides 20 5 --threads-per-worker 8  --n-parallel-bursts 2 --enable-gpu

vel=rasterio.open('timeseries/velocity.tif')
print(vel.shape)
plt.imshow(vel.read(1))

output: (928, 2241)

#Resample

os.makedirs('strided_geoms',exist_ok=True)
!/usr/bin/gdal_translate -outsize 2241 928 ../ISCE2-GPU_SLCstack/merged/geom_reference/lon.rdr.full.vrt strided_geoms/lon.strided.vrt
!/usr/bin/gdal_translate -outsize 2241 928 ../ISCE2-GPU_SLCstack/merged/geom_reference/lat.rdr.full.vrt strided_geoms/lat.strided.vrt

#Visualization

fig,axes=plt.subplots(ncols=2)
lon_strided=rasterio.open('strided_geoms/lon.strided.vrt')
lat_strided=rasterio.open('strided_geoms/lat.strided.vrt')
show(lon_strided,ax=axes[0])
show(lat_strided,ax=axes[1])

axes[0].set_title('Longitude')
axes[1].set_title('Latitude')

output:

#Utilize the warp function

def geocodeUsingGdalWarp(infile, latfile, lonfile, outfile,
                         insrs=4326, outsrs=None,
                         spacing=None, fmt='GTiff', bounds=None,
                         method='near'):
    '''
    Geocode a swath file using corresponding lat, lon files
    '''
    sourcexmltmpl = '''    <SimpleSource>
      <SourceFilename>{0}</SourceFilename>
      <SourceBand>{1}</SourceBand>
    </SimpleSource>'''
    
    driver = gdal.GetDriverByName('VRT')
    tempvrtname = 'geocode.vrt'
    inds = gdal.OpenShared(infile, gdal.GA_ReadOnly)
    
    tempds = driver.Create(tempvrtname, inds.RasterXSize, inds.RasterYSize, 0)
    
    for ii in range(inds.RasterCount):
        band = inds.GetRasterBand(1)
        tempds.AddBand(band.DataType)
        tempds.GetRasterBand(ii+1).SetMetadata({'source_0': sourcexmltmpl.format(infile, ii+1)}, 'vrt_sources')
  
    sref = osr.SpatialReference()
    sref.ImportFromEPSG(insrs)
    srswkt = sref.ExportToWkt()

    tempds.SetMetadata({'SRS' : srswkt,
                        'X_DATASET': lonfile,
                        'X_BAND' : '1',
                        'Y_DATASET': latfile,
                        'Y_BAND' : '1',
                        'PIXEL_OFFSET' : '0',
                        'LINE_OFFSET' : '0',
                        'PIXEL_STEP' : '1',
                        'LINE_STEP' : '1'}, 'GEOLOCATION')
    
    band = None
    tempds = None 
    inds = None
    
    if spacing is None:
        spacing = [None, None]
    
    warpOptions = gdal.WarpOptions(format=fmt,
                                   xRes=spacing[0], yRes=spacing[0],
                                   dstSRS=outsrs, outputBounds = bounds, 
                                   resampleAlg=method, geoloc=True)
    gdal.Warp(outfile, tempvrtname, options=warpOptions)

#Create and visualize geocoded output

geocodeUsingGdalWarp(infile='timeseries/velocity.tif',
                     latfile='strided_geoms/lat.strided.vrt',
                     lonfile='strided_geoms/lon.strided.vrt',
                     outfile='timeseries/velocity_geocoded.tif')

output:

[Shima-shoki]

@scottstanie
However, the InSAR seems executed not only over the land, but also over the ocean, so I'm little bit confused of interpreting the result 😅 Maybe I'd review the settings in dolphin, or in stackSentinel.py of the ISCE2...

[Shima-shoki]

@scottstanie
Oops, I should have used PS mask to get the reliable deformation result!
The outcome looks very cool!

vel=rasterio.open('timeseries/velocity.tif')
ps_pixels=rasterio.open('PS/ps_pixels_looked.tif').read(1)

vel_ps_masked=np.where(ps_pixels==1,vel.read(1),np.nan)
with rasterio.open('timeseries/velocity_ps_masked.tif','w',**vel.meta) as dst:
  dst.write(vel_ps_masked,1)

print(vel.shape)
plt.imshow(vel_ps_masked)

[Shima-shoki]

@scottstanie
A small question about the unit in the output. Based on the documentation, the unit of the final map might be rad/year. Therefore, if the "velocity" value at a certain place is -8, derived from Sentinel-1 (λ = 5.6 cm), the actual LOS change is -8/(2*pi) * λ = -7.1 cm/year. Is my interpretation correct...?

[scottstanie]

close! you need $4 \pi / \lambda$ , since the radar is two-way, so it will be half that velocity in cm. And you can double check the time series rasters with rio info --units ... . But I think you are correct