[Bug]: PRM.py

[kjz1997]

Describe the bug
hello professor,
when I use the read_SLC_int function of the PRM.py, it report the error, here is the error File "c:\Users\kjzha\.conda\envs\pygmt\Lib\runpy.py", line 198, in _run_module_as_main return _run_code(code, main_globals, None, ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "c:\Users\kjzha\.conda\envs\pygmt\Lib\runpy.py", line 88, in _run_code exec(code, run_globals) File "c:\Users\kjzha\.vscode\extensions\ms-python.debugpy-2024.8.0-win32-x64\bundled\libs\debugpy\adapter/../..\debugpy\launcher/../..\debugpy\__main__.py", line 39, in <module> cli.main() File "c:\Users\kjzha\.vscode\extensions\ms-python.debugpy-2024.8.0-win32-x64\bundled\libs\debugpy\adapter/../..\debugpy\launcher/../..\debugpy/..\debugpy\server\cli.py", line 430, in main run() File "c:\Users\kjzha\.vscode\extensions\ms-python.debugpy-2024.8.0-win32-x64\bundled\libs\debugpy\adapter/../..\debugpy\launcher/../..\debugpy/..\debugpy\server\cli.py", line 284, in run_file runpy.run_path(target, run_name="__main__") File "c:\Users\kjzha\.vscode\extensions\ms-python.debugpy-2024.8.0-win32-x64\bundled\libs\debugpy\_vendored\pydevd\_pydevd_bundle\pydevd_runpy.py", line 321, in run_path return _run_module_code(code, init_globals, run_name, ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "c:\Users\kjzha\.vscode\extensions\ms-python.debugpy-2024.8.0-win32-x64\bundled\libs\debugpy\_vendored\pydevd\_pydevd_bundle\pydevd_runpy.py", line 135, in _run_module_code _run_code(code, mod_globals, init_globals, File "c:\Users\kjzha\.vscode\extensions\ms-python.debugpy-2024.8.0-win32-x64\bundled\libs\debugpy\_vendored\pydevd\_pydevd_bundle\pydevd_runpy.py", line 124, in _run_code exec(code, run_globals) File "F:\下载文件\pygmtsar-pygmtsar2\pygmtsar-pygmtsar2\pygmtsar\temp2.py", line 44, in <module> block = dask.array.from_delayed( ^^^^^^^^^^^^^^^^^^^^^^^^ TypeError: from_delayed() missing 1 required positional argument: 'shape'
It seems that the return value of the read_SLC-block function does not meet the requirements of dask.array.for_delayed, so I change the code, just like values = dask.delayed(read_SLC_block(slc_file, start, stop))(2 * (stop - start)) block = dask.array.from_delayed( values, shape=(2 * (stop - start),), dtype=np.int16, ) then the code works

here is my all code :

import dask.delayed
import pygmt
import numpy as np
import xarray as xr
import numpy as np
import dask, dask.array
import os
import warnings

def read_SLC_block(slc_filename, start, stop):
return np.memmap(
slc_filename,
dtype=np.int16,
mode="r",
offset=start * 4,
shape=(2 * (stop - start),),
)
grd = r"F:\2018009_2018021\corr.grd"
slc_file = r"F:\raw\S1_20180110_ALL_F2.SLC"
scale=2.5e-07
xdim = 25272
ydim = 12192
chunksize = 2048
blocksize = chunksize * xdim
blocks = int(np.ceil(ydim * xdim / blocksize))
res = []
ims = []
for i in range(blocks):
start = i * blocksize
stop = min((i + 1) * blocksize, ydim * xdim)
values = dask.delayed(read_SLC_block(slc_file, start, stop))(2 * (stop - start))
block = dask.array.from_delayed(
values,
shape=(2 * (stop - start),),
dtype=np.int16,
)
res.append(block[::2])
ims.append(block[1::2])
del block
re = dask.array.concatenate(res).reshape((-1, xdim))[:ydim, :]
im = dask.array.concatenate(ims).reshape((-1, xdim))[:ydim, :]
del res, ims
coords = {"y": np.arange(ydim) + 0.5, "x": np.arange(xdim) + 0.5}
re = xr.DataArray(re, coords=coords).rename("re")
im = xr.DataArray(im, coords=coords).rename("im")
da = scale * (xr.merge([re, im]).astype(np.float32))

but I found the result is strange here is my result:

<xarray.Dataset> Size: 2GB
Dimensions: (y: 12192, x: 25272)
Coordinates:

• y (y) float64 98kB 0.5 1.5 2.5 3.5 ... 1.219e+04 1.219e+04 1.219e+04
• x (x) float64 202kB 0.5 1.5 2.5 3.5 ... 2.527e+04 2.527e+04 2.527e+04
  Data variables:
  re (y, x) float32 1GB dask.array<chunksize=(2048, 25272), meta=np.ndarray>
  im (y, x) float32 1GB dask.array<chunksize=(2048, 25272), meta=np.ndarray>
  it seems that the reuslt is not the SLC`s value, it's more like an offset of row and column numbers. can you give me some suggestion? look forward to your reply

[AlexeyPechnikov]

I’m not sure why the original code does not work for you, but the result looks okay. There are pixel-wise coordinates and the real and imaginary parts of the complex values stored in the SLC. The output is a lazy Xarray dask dataset, so you can get the values like da.re.values.

[kjz1997]

I’m not sure why the original code does not work for you, but the result looks okay. There are pixel-wise coordinates and the real and imaginary parts of the complex values stored in the SLC. The output is a lazy Xarray dask dataset, so you can get the values like da.re.values.

I use the da.re.values command, but it report the error, here is the error:
print(da.re.values) ^^^^^^^^^^^^ File "C:\Users\kjzha\.conda\envs\pygmt\Lib\site-packages\xarray\core\dataarray.py", line 786, in values return self.variable.values ^^^^^^^^^^^^^^^^^^^^ File "C:\Users\kjzha\.conda\envs\pygmt\Lib\site-packages\xarray\core\variable.py", line 540, in values return _as_array_or_item(self._data) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\kjzha\.conda\envs\pygmt\Lib\site-packages\xarray\core\variable.py", line 338, in _as_array_or_item data = np.asarray(data) ^^^^^^^^^^^^^^^^ File "C:\Users\kjzha\.conda\envs\pygmt\Lib\site-packages\dask\array\core.py", line 1693, in __array__ x = self.compute() ^^^^^^^^^^^^^^ File "C:\Users\kjzha\.conda\envs\pygmt\Lib\site-packages\dask\base.py", line 376, in compute (result,) = compute(self, traverse=False, **kwargs) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\kjzha\.conda\envs\pygmt\Lib\site-packages\dask\base.py", line 662, in compute results = schedule(dsk, keys, **kwargs) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\kjzha\.conda\envs\pygmt\Lib\site-packages\dask\array\chunk.py", line 420, in getitem result = obj[index] ~~~^^^^^^^ TypeError: tuple indices must be integers or slices, not tuple

[AlexeyPechnikov]

Try da[‘re’].values

[kjz1997]

Try da[‘re’].values

It reported the same error

[AlexeyPechnikov]

To verify your installation, run any of the provided PyGMTSAR projects. If everything functions correctly without any code modifications, your installation and function usage are properly set up. If issues arise, there may be a problem with your Python installation or function usage.

[kjz1997]

To verify your installation, run any of the provided PyGMTSAR projects. If everything functions correctly without any code modifications, your installation and function usage are properly set up. If issues arise, there may be a problem with your Python installation or function usage.

I have solved this bug, I changed the dask.array.from_delayed to dask.array.from_array, I compared the results of the two, they are the same, and the da.re.values was workable, thanks for your reply. But I still have some confusion, Why are both the azimuth and range coordinates offset by 0.5？

[AlexeyPechnikov]

The coordinate (0.5, 0.5) represents the center of a pixel with boundaries defined by the corners (0,0), (0,1), (1,1), and (1,0).

[kjz1997]

The coordinate (0.5, 0.5) represents the center of a pixel with boundaries defined by the corners (0,0), (0,1), (1,1), and (1,0).

ok, get it, thank you

[kjz1997]

I am studying the principle of InSAR by your code, yesterday I read the .SLC file under your guidiance, when I use the real value and image value to generate amplitude(sqrt(real2 + image2)), then show it, but the picture is too bad, here is the picture, can you give me suggestion? If you have offended me, please forgive me

here is the code:
amp_filt = np.sqrt((da.re.values**2 + da.im.values ** 2)) amp_filt[amp_filt == 0] = np.nan plt.imshow(amp_filt, cmap="gray", interpolation='nearest') plt.colorbar() plt.show()

[AlexeyPechnikov]

That’s ok for a linear scale because some pixels have much higher amplitudes. By the way, nodata pixels are not exactly zero.

[kjz1997]

Thank you,I solved it through histogram equalization

Have I noticed that the order of azimuth is reversed?

[AlexeyPechnikov]

Have I noticed that the order of azimuth is reversed?

How do you think what is the difference between ascending and descending orbits?:)

[kjz1997]

Have I noticed that the order of azimuth is reversed?

How do you think what is the difference between ascending and descending orbits?:)

When ascending, the radar image is inverted vertically, and when descending, the radar image is inverted horizontally. So when descending, do we need to change the order of distance? I am confused

[AlexeyPechnikov]

Of course, not. See the interactive PyGMTSAR notebooks for radar coordinates map examples. And you cannot change the radar coordinates because these are corresponding to the flight geometry.

[kjz1997]

ok，thanks for you, I will see the book

…

On 29 Jul 2024, at 00:45, Alexey Pechnikov ***@***.***> wrote:  Of course, not. See the interactive PyGMTSAR notebooks for radar coordinates map examples. And you cannot change the radar coordinates because these are corresponding to the flight geometry. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.

[AlexeyPechnikov]

Radar data is always processed on a uniform radar coordinate grid. While it is possible to find corresponding geographic coordinates for every radar pixel, these coordinates are uneven and cannot be processed directly. This means that only the processed results can be geocoded and exported.

[kjz1997]

In your book, I find that each pixel has its own incidence, so does each pixel also have its own spatial baseline?

[AlexeyPechnikov]

Baseline is satellite geometry characteristic so it is the same for all the pixels. Satellite look vector differ for every pixel and depending of the satellite position and topography.