AlexeyPechnikov / pygmtsar

PyGMTSAR (Python InSAR): Powerful and Accessible Satellite Interferometry
http://insar.dev/
BSD 3-Clause "New" or "Revised" License
432 stars 97 forks source link

[Feature]: OpenStreetMap Roads Mask #165

Open ClassiCcOwl opened 2 months ago

ClassiCcOwl commented 2 months ago

In S1A_Stack_CPGF_T173_roads notebook there is OpenStreetMap Roads , when in try to achive same result with PyGMTSAR 2 i get stuck in part where

sbas.unwrap_parallel(pairs, mask=osmmask)

can we have a updated version for that notebook with PYGMTSAR

AlexeyPechnikov commented 2 months ago

You’re trying to use outdated functions. I previously shared an example for sparse data processing, see the explanation at https://www.patreon.com/posts/pygmtsar-106488442 However, things are much simpler now.

ClassiCcOwl commented 2 months ago

I have phases and corrs gecoded using sbas.ra2ll() which is is provided in v2, and then i was able to get osmmask = sbas.ll2ra(osmmask_ll) ,i only need help on how to apply them together, in older version you have used masked paramater in unwrap_parallel function, but i couldn't find the equal function or parameter in new version.

AlexeyPechnikov commented 2 months ago

You can apply the mask in radar coordinates to your phase and correlation grids. After that, process them as usual.

ClassiCcOwl commented 2 months ago

You can apply the mask in radar coordinates to your phase and correlation grids. After that, process them as usual.

How can i apply that mask?

this is my corr and phase code. data = sbas.open_data() intensity = sbas.multilooking(np.square(np.abs(data)), wavelength=400, coarsen=(1,4)) phase = sbas.multilooking(sbas.phasediff(pairs), wavelength=400, coarsen=(1,4)) corr = sbas.correlation(phase, intensity) intf_filt = sbas.interferogram(sbas.goldstein(phase, corr, 32)) decimator = sbas.decimator() tqdm_dask(result := dask.persist(decimator(corr), decimator(intf_filt)), desc='Compute Phase and Correlation') corr60m, intf60m = result

tqdm_dask(phasefilts_ll := sbas.ra2ll(intf60m), desc='Geocoding') tqdm_dask(corrs_ll := sbas.ra2ll(corr60m), desc='Geocoding')

and this is my mask code

def dload_osm_roads(lon_min, lat_min, lon_max, lat_max):

OpenStreetMap graph calculations

import osmnx
# graph calculations
import networkx as nx

#north, south, east, west
G = osmnx.graph_from_bbox(lat_max, lat_min, lon_min, lon_max, network_type='drive')
osmroads = osmnx.graph_to_gdfs(G, nodes=False, edges=True)
return osmroads

workaround to run osmnx code in a separate process and fix the related Dask library error:

TypeError: _config_dns.._getaddrinfo() got an unexpected keyword argument 'family'

pathos isolate it well when joblib and multiprocessing do not

import pathos.multiprocessing as mp with mp.Pool(processes=1) as pool: osmroads = pool.apply(dload_osm_roads, args=(lon_min, lat_min, lon_max, lat_max))

add 90m buffer

osmroads_buffer = osmroads.to_crs("EPSG:32611").buffer(90).to_crs("EPSG:4326")

osmroads_buffer = osmroads.buffer(90)

osmmask_ll = sbas.as_geo(xr.ones_like(phasefilts_ll[0]).compute().rename({'lon': 'x', 'lat': 'y'})) osmmask_ll = osmmask_ll.rio.clip(osmroads_buffer).rename({'x': 'lon', 'y': 'lat'})

convert to radar coordinates

osmmask = sbas.ll2ra(osmmask_ll)

AlexeyPechnikov commented 2 months ago

Your osmmask should be the raster mask in radar coordinates. You might want to plot it to check. Simply multiply it by your interferogram or correlation raster, or use a condition like intf60m.where(np.isfinite(osmmask)) or intf60m.where(osmmask > 0) depending on the mask datatype (for float or boolean/binary masks).

ClassiCcOwl commented 2 months ago

when i use

sbas.plot_interferograms(intf60m.where(osmmask > 0), cols=3, size=3, caption='Phase, [rad]') plt.savefig('Phase masked, [rad].jpg')

it never stops calculating

can you update the notebook with new version please?

AlexeyPechnikov commented 2 months ago

If your osmmask is a lazy object, materialize it with osmmask.compute() or by using persist() with a progress bar, or by using PyGMTSAR’s sync_* functions. Additionally, you can use a materialized intf60m if you still encounter performance issues.

ClassiCcOwl commented 2 months ago

Thanks i will try them

ClassiCcOwl commented 2 months ago

image

still wasn't able to fix it

AlexeyPechnikov commented 2 months ago

You don’t need to focus on ‘INFO’ messages, as these are typically logged for the benefit of Dask developers and are meant for informational purposes only. If seeing these messages is bothersome, you can opt to run the Dask cluster from the command line outside of the notebook and connect to it from within your notebook. This approach helps prevent the display of such messages.

ClassiCcOwl commented 2 months ago

When i use phase.where(osmask) i get no result in plots because it says there is no numeric to plot, apparently the osmmask that i have and phases do not fit or work with each other like they were used to. What alternative method do we have for sbas.unwrap_parallel(pairs, mask=osmmask) in new version of pygmtsar, sir?

AlexeyPechnikov commented 2 months ago

As I mentioned earlier, you need to verify whether your mask is in float or boolean format and apply it accordingly. Please refer to the examples provided above for guidance.