Otsus method

I reached the point has mentioned via this post in the last few days,

I found that unfortunately there’s still no method via SNAP, to calculate the intensity threshold of VH, VV,

That’s why I edited and created the Otsu’s method via skimage library to be compatible with GRDHR images, and full preparing the intensity VH, VV, and also getting the threshold is available within my github repository.

2 Likes

sounds cool - I would love to see an example for its use.

Dear Andy,
An example is already mentioned via the repository,

dirpath = ‘D://WATER_MAP_TEST/’
vv = dirpath + VV/S1A_IW_GRDH_1SDV_20210324T191554_VV_NR_Orb_TC_res.tif’
vh = dirpath + ‘VH/S1A_IW_GRDH_1SDV_20210324T191554_VH_NR_Orb_TC_res.tif’
Otsus_VH_VV_thr(vv)
layer_count : 1 layer_shape : (21494, 25689)
layer_type : (‘float32’,)
Lyer_NAN : (None,)
Layer_shape_after_converted_to_unit8 (1, 21494, 25689)
The layer threshold is : 21
Total time is: 00:07:01.82

image

sorry, I meant the result.

Yes, you’re right I miss-interpreted your previous comments, the results of thr. are required in order to create water extent mapping as it’s explained here, make_water_map

This is an example,

import asf_tools
from asf_tools.water_map import make_water_map

dirpath = ‘D:/WATER_MAP_TEST/’

vh = dirpath + ‘VH/S1A_IW_GRDH_1SDV_20210324T191554_VH_NR_Orb_TC_res_Ltodb.tif’

vv = dirpath + ‘VV/S1A_IW_GRDH_1SDV_20210324T191554_VV_NR_Orb_TC_res_Ltodb.tif’

out_raster = dirpath + ‘water_map_20210324T191554.tif’

make_water_map(
out_raster,
vv,
vh,
hand_raster=None,
tile_shape=(200, 200),
max_vv_threshold=21,
max_vh_threshold=18,
hand_threshold=15.0,
hand_fraction=0.8,
membership_threshold=0.45
)

I hope now the foregoing answers your comment!

1 Like

I faced similar kind of issue last time, I am still searching for some proper solution.

Hi Orland,

Seemingly you just joined the step_forum, I’d like to welcome you in the world of Space and Remote sensing,

Would you please to clarify which kind of problem do you face up?, your comment is vague to me!

Dear @falahfakhri ,
thank you to share this code. Indeed it is very useful.
I am working with the same kind of data (IW GRDH) and I would like to apply the Otsu filter to separate land from the sea and extract a coastline.
I have a question about your code.
Firstly, is there any reason in your preprocessing any calibrations is applied? It seems you are working just with the DN and then after the Terrain correction you transform your data from linear to db space.
Secondly, before applying the conversion with

img = img_as_ubyte(vvvh)

, you divide your image by 255. Are you doing that to scale it? If yes, shouldn’t it be something like

img= ((img- np.min(img)) * (1/(np.max(img) - np.min(img)) * 255)).astype(‘uint8’)

I’m asking because my data after my preprocessing (that include several steps like calibration, Terrain Flattening, etc.) has values in the range in dB [-54.40889358520508, 16.46017837524414].
When I divide them by 255, they are in the range [-0.21336821013805912, 0.06454971911860448]
and so when I transform the data with img_as_ubyte they become in the range [0,1] and the filter doesn’t work.
On the other side, if I scale using the formula above, I will get values between [0, 255] and in this case, the filter gives a reasonable output.
I would like to have your thoughts about this. Any guidance is warmly appreciated.

S Savastano

1 Like

Hi and thanks a lot for your comments and questions.

First of all I already shared the graph you can update it and follow the same steps. The last step is to follow what I mentioned.

– Thresholding works best when backscatter tiles are provided on a decibel scale to get Gaussian distribution that is scaled to a range of 0-255, and performed on a small tile that is likely to have a transition between water and not water.

– In SNAP after the graph has been run and the intesity of the two pols. has been prepared, in the product explorer expand the bands, right click and select Linear/to from dB. Extract the new dB band and from File tab export it as BigGeotiff.

No.

Hi @falahfakhri ,
so do you mean your data after the conversion in dB is already scaled in the range [0, 255]?

Once the processed tile converted to decibel the backscatter are then scaled to get Gaussian distribution that is scaled to a range of 0-255.

Hope it’s clear now.

Hi,
it’s clear what you mean but I don’t understand when your data is scaled. After running your graph, you get a .dim file with 2 bands. Then you convert in SNAP from linear to dB and finally, you export it as a Geotiff image.
Now my question is: if your data is in dB, how you can be sure they are in the range [0,255]?
is it the rasterio command in python that scales your dB values in that range?
Sorry to bother you, but I would like to use your function but I would like to understand how to suit it to my data.

No worries,

Just check your result via SNAP pixel value to be sure of your result.

this is my data after conversion in dB

So, when and how can I scale it to be sure to be in the range [0, 255] and apply Otsu in the best input conditions?

Right click on the product and then convert to db. Unfortunately I still don’t have SNAP in this computer to show you by screenshot.

Hi @falahfakhri ,
maybe there is a misunderstanding.
My data after the preprocessing is already in dB as you can see from the plot above.
I am asking, when and how do you scale your data in dB to have them in the range [0, 255] that is a compulsory requirement to apply Otsu?
I don’t see this scaling once you applied your preprocessing.

Thank you

Hi,
if you give a quick glance at the code you could find the the threshold_otsu modules is already imported please see the screenshot below

image

Then within the code the local and the global otsu are applied, please see the screenshot below.

image

You don’t need to scale your data.

I hope I successfully covered your requirement.