Help starting with PyroSAR

[SantaTitular]

Hi,

I'm just starting to use pyroSAR after searching for a simpler way to process Sentinel S1 images in python. I started by just doing some standard processing but I'm not sure how to include the geometry (wkt) and the projection (proj) as if using snappy. I tried looking through the function itself but the t_srslink is broken unfortunately. The code below naturally outputs an error: OSError: file does not exist in the sub = sub_parametrize(scene=id, geometry=shapefile, offset=offset, buffer=0.01), shp = Vector(geometry).

wkt = 'POLYGON((19.295571567283318 43.75986021410339, 19.293968420950073 43.759242194967975, \
19.288328627371172 43.7581532581082, 19.286606237262227 43.75785887570865, \
19.28440871687089 43.75678787464691, 19.282673948937227 43.756404046881535, \
19.280178700089042 43.75560522857542, 19.277021698089534 43.7533436552217, \
19.275516711676158 43.75171595614192, 19.275648134538415 43.75041233793973, \
19.274814221707484 43.74898099255665, 19.273359746684783 43.74787109709873, \
19.27198229230851 43.74725418142482, 19.26825816651079 43.74723157460489, \
19.262711461336814 43.74803405929451, 19.25806578020773 43.749332419458206, \
19.251175882716723 43.752117208490326, 19.2519962717524 43.75297927471196, \
19.25532370925354 43.75181462755526, 19.25906504891751 43.75045575733781, \
19.262007119781806 43.74976883625828, 19.263486703638296 43.74967431266951, \
19.26663866444263 43.748380844055504, 19.26881463420939 43.748253842277116, \
19.271066206745964 43.74841579293915, 19.272269863969786 43.74909182415384, \
19.273283239646606 43.750304778053945, 19.27502656136605 43.75248081272469, \
19.276712077409602 43.75463398741618, 19.279689628172484 43.75670609886505, \
19.281267857043577 43.7572830742399, 19.282084907473866 43.758100661421324, \
19.28418263838442 43.75911170596844, 19.285504881814557 43.75982192677596, \
19.287884999799278 43.7608787870774, 19.290834466804345 43.76076046044949, \
19.292899545616095 43.761095905420454, 19.294267802847358 43.76147076065607, \
19.294267802847358 43.76147076065607, 19.295571567283318 43.75986021410339))'
## UTM projection parameters
proj = '''PROJCS["UTM Zone 4 / World Geodetic System 1984",GEOGCS["World Geodetic System 1984",DATUM["World Geodetic System 1984",SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],UNIT["degree", 0.017453292519943295],AXIS["Geodetic longitude", EAST],AXIS["Geodetic latitude", NORTH]],PROJECTION["Transverse_Mercator"],PARAMETER["central_meridian", -159.0],PARAMETER["latitude_of_origin", 0.0],PARAMETER["scale_factor", 0.9996],PARAMETER["false_easting", 500000.0],PARAMETER["false_northing", 0.0],UNIT["m", 1.0],AXIS["Easting", EAST],AXIS["Northing", NORTH]]'''
input_S1_files = sorted(list(iglob(join(path, '**', '*S1*.zip'), recursive=True)))

for i in input_S1_files:
    print(input_S1_files)
    geocode(i,
        outdir=outpath, test=True, t_srs=proj,
        polarizations=['VV','VH'],
        speckleFilter='Lee',
        removeS1ThermalNoise=True,
        allow_RES_OSV=True,
        shapefile=wkt)
    break

Thanks in advance! Tomás

[johntruckenbrodt]

Hi @SantaTitular I agree the documentation of the parameter shapefile is not very clear and the name itself is misleading. shapefile is a file format and the function accepts more than that. The parameter geometry of function snap.auxil.sub_parametrize (which is called by geocode) gives a better description. Furthermore, I see that the link to osgeo.osr.SpatialReference is broken, I guess this is the one you mean. I will update the documentation of the shapefile parameter accordingly.
Here's what you can do (using function spatialist.vector.wkt2vector):

from spatialist.vector import wkt2vector

for i in input_S1_files:
    print(i)
    with wkt2vector(wkt=wkt, srs=proj) as vec:
        geocode(i,
            outdir=outpath, test=True, t_srs=proj,
            polarizations=['VV','VH'],
            speckleFilter='Lee',
            removeS1ThermalNoise=True,
            allow_RES_OSV=True,
            shapefile=vec)    

Works for you?

[SantaTitular]

Hi @johntruckenbrodt

Thanks so much for the response! Yes, it was the link for osgeo.osr.SpatialReference that was broken.

You're right, after looking through the documentation of the sub_parametrize function, it became quite clear (great tip on the use of the wkt2vector function!). Now when I run the code I get an error from the intersection of scene and shapefile (image below). Maybe its an error from the very specific polygon I used or a difference in the Spatial reference (This is the image I tested in SNAP: S1A_IW_GRDH_1SDV_20221004T164201_20221004T164226_045295_056A44_9032).

Thanks again!

Quick edit: It does work without the subset (atleast its running and created the .xml file previously), so it is just this subset part.

[johntruckenbrodt]

Your AOI is located in the central Pacific near Samoa, the scene covers Bosnia Herzegovina. Looks like you need to adjust your WKT string. the coordinate convention is always X-Y(-Z) whereas you seem to have Y-X.

- POLYGON((19.295571567283318 43.75986021410339, 
+ POLYGON((43.75986021410339 19.295571567283318, 

[SantaTitular]

Thanks again for the help. Although I'm confused since I took the polygon's WKT from SNAP. Still, if I look up this polygon POLYGON ((19.279547 43.774936, 19.279547 43.794084, 19.304609 43.794084, 19.304609 43.774936, 19.279547 43.774936)) on a WKT visualization I get the coordinates I wanted.

I did try changing the X-Y on the code but ended with the same error!

[johntruckenbrodt]

Ah sorry. Swapping the coordinates would obviously get you to the Arabian peninsula and not to Samoa. There must be a different reason. Until I get this fixed I suggest you find a different way to create some vector file from your WKT string and use this file for processing.

[SantaTitular]

Hi again John,

I tried processing without the subset however, I ended with an error file with the following line on it: 'utf-8' codec can't decode byte 0xe3 in position 2027: invalid continuation byte. Maybe it is something I inputted wrong, but I just want to do a standard preprocessing with the set of images.

Edit: File name: S1A__IW___A_20221004T164201_Cal_NR_Orb_ML_TF_Spk_TC_dB_error

[johntruckenbrodt]

Okay I figured out the issue with the geometry. Your WKT string coordinates do not match the coordinate system of proj.
With this is works:

proj = 'EPSG:4326'

There is still a shift of the geometry relative to the river though:

[johntruckenbrodt]

Your last message looks a lot like this: https://github.com/johntruckenbrodt/pyroSAR/issues/265
Do you perhaps have any special characters in the file names written to your workflow?

[SantaTitular]

Hi John,

Thanks again for the help! Thats weird, If I remember correctly, I took proj's WKT from SNAP (Raster->Geometric->Reprojecting) but maybe before doing any shift in the image or drawing the polygon. I'll just use the projection from WKT online viz for future cases then.

Regarding the error message, contrary to #265 this error would appear at the end, at the write module. I imported the graph to SNAP as you suggested in #265 (although he did have '(' special characters in the path name) and changed directory and the output format:

BEAM-DIMAP

GeoTIFF

ENVI

However, if I incorporate the subset in the workflow (using the original, possible wrong path) I do get an output (albeit still not the one I wanted 😅)! So I guess it is more related to size of the image than anything else.

[johntruckenbrodt]

Mh strange. I think you're right, this looks a lot like you don't have enough space.
I am not sure I can follow, how is the output not what you wanted?

[SantaTitular]

Yeah, I guess 500GBs not enough to process a full image. Still, I can process with the subset so I'm ok with it. Since I used lineartodB the output image was pretty much black, whereas I just wanted the linear form like this:

I don't have much experience with this type of data but I think it look ok no? I'm still exploring more of the package and different processing steps and it is really helpfull to check the intermediate steps in SNAP. Thanks for the help again!

[johntruckenbrodt]

500 GB should be plenty to process the images so the error could be coming from something else.
I think you are interchanging linear and logarithmic (dB) scale. In log scale the image should be much better to view than in linear scale. Your image looks alright but without any scale it is hard to tell what you show.

[SantaTitular]

Hi John,

I'm back to using the pyroSAR library after taking time off but I can't seem to get the pyroSAR library working as previously. I lost the conda environment that I was using before and I can't seem to make a working one.

I started by creating an environment with python 3.6 and installing pyroSAR (using conda install --channel conda-forge pyrosar) but the code only supported python > version 3.10? But using an environment with version 3.10 or 3.8 I can't seem to configure snappy nor the correct dependencies (for instance, I get the message 'no module osgeo') Don't suppose you have a .yaml file to help me out?

Best, Tomás

**EDIT: I Installed SNAP v10.0.0 and now I'm getting the following error: AttributeError: SNAP could not be identified. If you have installed it please add the path to the SNAP executables (bin subdirectory) to the PATH environment. E.g. in the Linux .bashrc file add the following line: export PATH=$PATH:path/to/snap/bin"

**EDIT2: Persistence pays off, after updating everything and restarting the computer I managed to run the code that you suggestted above (below). However, when I put test = False to run the .xml workflow, I still have memory issues as before. Is there something I can change just to obtain the correct final image?

from spatialist.vector import wkt2vector

for i in input_S1_files:
    print(i)
    with wkt2vector(wkt=wkt, srs=proj) as vec:
        geocode(i,
            outdir=outpath, test=True, t_srs=proj,
            polarizations=['VV','VH'],
            speckleFilter='Lee',
            removeS1ThermalNoise=True,
            allow_RES_OSV=True,
            shapefile=vec)    

[johntruckenbrodt]

@SantaTitular great you were able to get it to run. Sorry for all the issues you're encountering.
The function snap.geocode has a parameter gpt_args, which lets you configure the resources used by SNAP. For example, you can do the following to adjust memory (32 GB), cache (22 GB) and number of threads (4): gpt_args = ["-J-Xmx32G", "-c", "22G", "-q", "4"]

[SantaTitular]

@johntruckenbrodt thanks for the tip. Yesterday I had it working perfectly (took a while to process all the images I had), but for some reason today I had the same error:

version = snap.get_version('microwavetbx')['version']

  File ~\anaconda3\envs\snap310\lib\site-packages\pyroSAR\examine.py:332 in get_version
    raise RuntimeError('cannot read version information from {}.\nPlease restart SNAP.'.format(fname))

RuntimeError: cannot read version information from C:\Users\tomas\AppData\Roaming\SNAP\var\log\messages.log.
Please restart SNAP.

I went to the code part raising the error and found the following:

I'm not sure what originated this problem, all I did different was open SNAP and check for updates (maybe installed something bad?). If you're available, we could discuss further via email or a quick chat as I'm really keen on using pyroSAR for my project!😁

*EDIT: Yeah, after I removed and renistalled SNAP, it went back to working. Maybe it had something to do with the plugin installer... Still, if you're available, I would really enjoy discussing some questions regarding exclusions layers, etc.

[johntruckenbrodt]

I have written a mail. Has it reached you?