Extract_information

[giannis]

Is it possible to use extract Pixel values but instead of writing them in txt store them in a variable e.g. list, array in python? My point is that i have multiple products and relevant coordinates (in some products i have two or more coordinates to extract) but i do not want to create many txt/.

[M.Zohary]

take a look
this may help
def product2array(src_product, band_names=[] ):
width = src_product.getSceneRasterWidth()
height = src_product.getSceneRasterHeight()
name = src_product.getName()
desc = src_product.getDescription()
if not band_names:
band_names = src_product.getBandNames()
# band_names = [b for b in band_names if “Sigma” in b]
print(‘We are in product2array extracting this {} \n from this {}’. format(band_names,name ))
print('product2array width:{} height:{} name:{} /n desc:{} /n band_names: {} '.format(width, height, name, desc, band_names))

‘’’ creating empty aprray ‘’’
filename = path.join(mkdtemp(), ‘dataStorageFile.dat’)
# Create a memmap with dtype and shape that matches our data:
data = np.memmap(filename, shape=(height, width, len(band_names)), dtype=np.float32, mode=‘w+’)
tempname = path.join(mkdtemp(), ‘tempStorageFile.dat’)
temp = np.memmap(tempname, shape=(height, width), dtype=np.float32, mode=‘w+’)
for b in range(len(band_names)):
# print(‘reading band:’, band_name)
print(’ reading: ', band_names[b])
band = src_product.getBand(band_names[b])
band.readPixels(0, 0, width, height, temp)
data[:, :, b] = temp
return data

[marpet]

You can use the PixEx operator, it writes to a text file. But you can easily read it in again.
Another option is you do the extraction your self. You can iterate over all bands of a product and get the value at the specific location. To convert your geo-location into a pixel location and read the value you can do:

gc = product.getSceneGeoCoding()
pixelPos = gc.getPixelPos(GeoPos(lat, lon), None)
value = band.getPixelDouble(pixel.Pos.x, pixelPos.y)

This is roughly what you need to do.

[giannis]

Trying to adopt your code

product=ProductIO.readProduct(folder+‘S3A_OL_2_WFR____20160430T093820_20160430T094020_20171030T225518_0119_003_307______MR1_R_NT_002.SEN3/xfdumanifest.xml’)
gc = product.getSceneGeoCoding()
pixelPos = gc.getPixelPos(GeoPos(40.36734009, 10.64342976), None)
band = product.getBand(‘CHL_OC4ME’)
value = band.getPixelDouble(pixelPos.x, pixelPos.y)
print value

but i get the following error File “E:/seo_dwarf/chl_scripts/extract_multiple_products1.py”, line 93, in
value = band.getPixelDouble(pixelPos.x, pixelPos.y)
RuntimeError: no matching Java method overloads found

Any ideas?

[marpet]

The variables pixelPos.x and pixelPos.yare floating point. And getPixelDouble() expects integers.
The values need to be cast to integer. This should do the trick.

value = band.getPixelDouble(int(pixelPos.x), int(pixelPos.y))

[giannis]

Thanks for your answers but now i get this error

value = band.getPixelDouble(int(pixelPos.x), int(pixelPos.y))
RuntimeError: java.lang.NullPointerException

Any ideas?

[marpet]

Oh, yes. My mistake. The getPixelDouble method expects that the data is already loaded.

Either do

band.loadRasterData() 

before calling getPixelDouble (once per band is sufficient) or

data = array.array('d', [0]) 
# or with numpy
# import numpy as np
# data =np.zeros(1)
readPixels(int(pixelPos.x), int(pixelPos.y), 1, 1, data)

[giannis]

Is there any possibility to use the float and not convert the coordinates to integers?

[marpet]

no, and it makes no difference.
You will get the same pixel value either if you use

x = 34.567
y = 123.456

or
x = 34
y = 123

[giannis]

Did what you said but i get the following error

Traceback (most recent call last):
File “<pyshell#3>”, line 1, in
value = band.getPixelDouble(int(pixelPos.x), int(pixelPos.y))
RuntimeError: java.lang.ArrayIndexOutOfBoundsException: 6027866

Any ideas?

[marpet]

Is the pixel not anymore in the bounds of the image?

catch the exception and check at which x,y coordinate this happens. Check also the image bounds. Similar to this.

try:
    value = band.getPixelDouble(int(pixelPos.x), int(pixelPos.y))
except RuntimeError:
    print "x-coord:", int(pixelPos.x)
    print "y-coord:", int(pixelPos.y)
    print "band width:", band.getRasterWidth()
    print "band height:", band.getRasterHeight()

[Pang]

Dear Marco,
There is an error from:

pixelPos = gc.getPixelPos(GeoPos(30.7633,90.7628), None)
RuntimeError: no matching Java method overloads found

These are whole code:

from snappy import ProductIO
product=ProductIO.readProduct(‘I:/Sentinel3/PixEx-input/gaican-0726.nc’)
gc = product.getSceneGeoCoding()
pixelPos = gc.getPixelPos(GeoPos(30.7633,90.7628), None)
band = product.getBand(‘rhwon_3’)
value = band.readPixels(pixelPos.x, pixelPos.y, 1, 1).getFloat(0)
print(value)

Thanks so much.

[marpet]

@Pang

pixelPos = gc.getPixelPos(GeoPos(30.7633,90.7628), None)

Can you split the line into two.
Then it will be clear if getPixelPos is the problem or the GeoPos creation.
Instead of None you can try an empty PixelPos().

For others this type of call seems to work:

Get a tile's geocoding data

I see, thanks for the answer. Yes, only one source. I was under the impression that the result given by getPixelPos() would be oriented relative to the full raster. For example, if a given raster is 1000x1000, and a given lat-long corresponds to: self.source_product.getSceneGeoCoding().getPixelPos(pos, None).getX() = 995 self.source_product.getSceneGeoCoding().getPixelPos(pos, None).getY() = 995 So at an index of (995, 995) in the original product’s raster. But if we use this index in a 100…

[Pang]

Dear Marco,
Thanks so much for your help.
Those codes can output data that I wanted (rhown_3).

from snappy import ProductIO
from snappy import PixelPos, GeoPos
import array
product=ProductIO.readProduct(‘I:/Sentinel3/PixEx-input/gaican-0726.nc’)
gc = product.getSceneGeoCoding()
pixelPos = gc.getPixelPos(GeoPos(30.7633,90.7628), None)
geo_pos = GeoPos(30.7633, 90.7628)
pixel_pos = gc.getPixelPos(geo_pos, PixelPos(0,0))
band = product.getBand(‘rhown_3’)
band.loadRasterData()
data = array.array(‘d’, [0])
value = band.readPixels(int(pixelPos.x), int(pixelPos.y), 1, 1, data)
print(data[0])

Then I want to batch them and obtain rhown1-21.

[Pang]

Dear Marco,
I want to know if I should modify the parameters?
My target is the lake of the Tibet Plateau, they are very clean. TSM and Chla concentrations are very low. There are high altitude and low atmospheric pressure. I have just tried to modify only Salinity is 1.1, Temperature is 10.1, rhown_3 was 0.02956 (before) and 0.02926 (after). They are not very different, so if not necessary to modify the parameters?

[Pang]

Dear Marco,
In addition to C2RCC, what preprocessing is needed for OLCI L1 like radiative calibration? C2RCC seems to include all operations like radiation correction?
And Rrs=rhown/π right？How about rrs？Rrs=rrs？

Thanks a lot.

[marpet]

In some cases, it helps to apply vicarious calibration. This is explained an FAQ entry.
S3TBX FAQs - SNAP - Confluence (atlassian.net)
Otherwise, no pre-processing is required.

The π is considered by SNAP.

[Pang]

Dear Marco,
Thanks for your reply.
But my rhown and rrs were not from SNAP, they were processed by C2RCC in python. So I should used rrs rather rhown right?

And what is your idea about this question?

Dear Marco,
I want to know if I should modify the parameters?
My target is the lake of the Tibet Plateau, they are very clean. TSM and Chla concentrations are very low. There are high altitude and low atmospheric pressure. I have just tried to modify only Salinity is 1.1, Temperature is 10.1, rhown_3 was 0.02956 (before) and 0.02926 (after). They are not very different, so if not necessary to modify the parameters?

Thanks a lot.

[marpet]

I’m not a remote sensing scientist. Just a software developer. So, I fear this is not my domain and I shouldn’t reply.

Maybe @abruescas can provide more information.

[abruescas]

I think you can use the Rrs (sr-1), but we aware of the units you use before comparing with in situ data.
In any case, if your lake is very clear, you can also check the ESA L2 standard product directly. C2RCC is designed for complex waters, and maybe you will not obtain good results.