Get latitude/longitude values in python operator

development python
1

Hello,

I am currently programming a python operator for snap.
In the method “computeTileStack” I get the data from the bands in tiles (getSourceTile) and read the samples (getSamplesFloat) from there. I convert this to a numpy array for convenience. This is working fine.

Additionally, I need the latitude and longitude values for each pixel (needed for lookup in a lookup table). So I would like to have a numpy array of the same size as the samples array with the latitute data and one with the longitute data.

My first try was to use a bandmath expression to add a new band to the source product that contains this information (expression = ‘LAT’ or ‘LON’). Then I also can use getSourceTile and so on. But these additional bands should not be added to the source or target product. Is there another way to get this information without creating new bands?

Thanks for any tipps.
Kind regards
Helga

2

Hi Helga,

actually the official intended way to get the lat/lon values is via the Geocoding.

from snappy import (GeoPos, PixelPos)
gc = sourceProduct.getSceneGeoCoding()
#or in multi-resolution case use the geo-coding of the raster
gc = sourceProduct.getRasterDataNode(“name”).getGeoCoding()
geoPos = gc.getGeoPos(PixelPos(x, y), None)
lat = geoPos.getLat()
lon = geoPos.getLon()

You would iterate over your pixel indices (x,y) and fill an array. But this can be tedious and actually you already found the solution. At least halfway. You can use setOwner(). This way the new band is owned by the source product but is not part of it.

virtLat = VirtualBand(“__lat”, ProductData.TYPE_FLOAT32, sourceProduct.getSceneRasterWidth(), sourceProduct.getSceneRasterHeight(), “LAT”)
virtLat.setOwner(sourceProduct)
virtLon = VirtualBand(“__lon”, ProductData.TYPE_FLOAT32, sourceProduct.getSceneRasterWidth(), sourceProduct.getSceneRasterHeight(), “LON”)
virtLat.setOwner(sourceProduct)

3 Likes
3

Thanks. This absolutly solved the problem.

4

Hi there,
I’m using bandmath expression to get the information about LAT and LON and add them as bands in a product.
My problem is that the operator seems to work well until a point where (don’t know for which reason) the pixels get assigned to the same value.
Here is my code:

BandDescriptor = jpy.get_type('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor')
        
lat = BandDescriptor()
lat.name = 'Latitude'
lat.type = 'float64'
lat.expression = 'LAT'
newBand = jpy.array('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor', 1)
newBand[0] = lat
parameters = HashMap()
parameters.put('targetBands', newBands)
prod = GPF.createProduct('BandMaths', parameters, source)

If I plot the pixel values of the ‘Latitude’ band inside prod I get the following results:

Cattura.PNG1235×803 18.7 KB

As can be seen, in the bottom area the pixels are all of the same value (41.31)

If I try doing the process on the same image through Band Maths operator in SNAP GUI, I get the correct output:

Cattura1.PNG1100×792 2.41 KB

What is wrong with my BandMaths procedure for getting LAT value? Is there a wat to get the correct values through snappy?
Thank you!

Kind regards,
Davide

5

On which data are you letting this run?

6

The code snippet I wrote above is the ADD LAT and LON BANDS TO THE PRODUCTS part of the stacking code I created (You can find it here: Error on writing target product).
In that post, that part is commented because the behaviour I got from the bandmaths operator wasn’t as expected (as I was explaining in my post above).

7

It seems you dug out a bug which was undiscovered for years.
I’ve created an entry for this in our issue tracker.

You can workaround if you use the names of the tie-point grids directly, I guess you still have them in your sources.
Just replace in the expression ‘LON’ by ‘longitude’ and ‘LAT’ by ‘latitude’.

Thanks for your report and I hope the workaround will work.

8

Thank you for your reply Marco, I understand.
I managed to solve my problem without having the Lat and Lon bands actually so I don’t need this anymore, but I’ve tried your workaround and unfortunately it doensn’t work :confused: it raises the error: RuntimeError: org.esa.snap.core.gpf.OperatorException: Could not parse expression: 'latitude'. Undefined symbol 'latitude'

9

Just wanted to note here that (x, y) in PixelPos(x, y) should be interpreted as positions along the pixel matrix in units of of pixel width and height, respectively.

For instance, gc.getGeoPos(PixelPos(0, 0), None).getLon() would return the longitude at the upper left corner of the pixel matrix, which is also the upper left corner of the pixel of indexes (0, 0). This is not the longitude at the centre of such pixel! If it is the longitude of such centre what is wanted, one would have to use gc.getGeoPos(PixelPos(0.5, 0.5), None).getLon() instead.

And, for instance, if it is the longitude at the centre of the pixel at the lower right corner of the pixel matrix what is wanted, then gc.getGeoPos(PixelPos(width - 0.5, height - 0.5), None).getLon() would need to be used. Here, width and height are the width (number of columns) and height (number of rows) of the pixel matrix, respectively.