Search code examples
pythongeolocationsentinel1snap-pythonsentinel2

Converting geographical coordinates (long, lat) to pixel position (x,y) for a Sentinel-1 SAR image


How can i get the (x, y) pixel position from the geographical coordinates in a Sentinel-1 Synthetic Aperture Radar (SAR) satellite image?

I can access a downloaded image info sg as

from snappy import ProductIO

path='path_name'
product = ProductIO.readProduct(path)
sg = product.getSceneGeoCoding()

But how can I get (x, y) pixel position for a desired latitude and longitude using ESA's snap engine within Python?


Solution

  • Using the custom function below, we can easily convert any (latitude, longitude) to it' s (x, y) position in image, provided that latitude and longitude is within the limits of our product.

    from snappy import GeoPos
    def XY_from_LatLon(ProductSceneGeoCoding, latitude, longitude):
        #From Latitude, Longitude satellite image (SAR), get the x, y position in image
        pixelPos = ProductSceneGeoCoding.getPixelPos(GeoPos(latitude, longitude), None)
        x = pixelPos.getX()
        y = pixelPos.getY()
        if str(x)=='nan':
            raise ValueError('Latitude or Longitude out of this product')
        else:
            return x, y
    

    UPD: Updated function below should work with more snap versions

    import jpy
    import snappy
    def XY_from_LatLon(ProductSceneGeoCoding, latitude, longitude):
        geoPosType = jpy.get_type('org.esa.snap.core.datamodel.PixelPos')
        geocoding = ProductSceneGeoCoding.getSceneGeoCoding()
        pixel_pos = geocoding.getPixelPos(snappy.GeoPos(latitude, longitude), geoPosType())
        if str(pixel_pos.x)=='nan':
            raise ValueError('Latitude or Longitude out of this product')
        else:
            return int(np.round(pixel_pos.x)), int(np.round(pixel_pos.y))
    

    E.g. for the given product bellow (as we can see in scihub it's a product in southern Greece), we can get (x, y) position in image of Athens coordinates (latitude=37.9838, longitude=23.7275) as

    productname : S1A_IW_GRDH_1SDV_20170821T162310_20170821T162335_018024_01E414_C88B

    path='path to S1A_IW_GRDH_1SDV_20170821T162310_20170821T162335_018024_01E414_C88B.SAFE'
    product = ProductIO.readProduct(path)
    sg = product.getSceneGeoCoding()
    x, y = XY_from_LatLon(sg, 37.9838, 23.7275)
    x, y
    # (13705.242822312131, 14957.933651457932)