Get output coordinates using sentinelhubpy

Hi everyone,

I’m using the Python API of Sentinelhub for a couple of weeks and made my first requests using a handmade function, here is an example:

from sentinelhub import SHConfig
from sentinelhub import MimeType, CRS, BBox, SentinelHubRequest, SentinelHubDownloadClient, \
     bbox_to_dimensions, DownloadRequest, DataCollection
import requests


def link_to_text(link):
    f = requests.get(link)
    return f.text

#NDVI
escript_NDVI=link_to_text("https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/ndvi/script.js")

CLIENT_ID = my_id
CLIENT_SECRET = my_secret

config = SHConfig()

resolution=40

if CLIENT_ID and CLIENT_SECRET:
    config.sh_client_id = CLIENT_ID
    config.sh_client_secret = CLIENT_SECRET

if config.sh_client_id == '' or config.sh_client_secret == '':
    print("Warning! To use Sentinel Hub services, please provide the credentials (client ID and client secret).")

data_folder=my_folder
time_interval=('2020-07-01', '2020-07-02')

coords = [122.638, 64.229, 124.262, 65.043]

def sentinelhub_request(data_folder=data_folder,time_interval=time_interval,coords=coords,
                        evalscript=escript_active_fire,
                        save_data=False):
    
    loc_bbox = BBox(bbox=coords, crs=CRS.WGS84)
    loc_size = bbox_to_dimensions(loc_bbox, resolution=resolution)
    
    request_all_bands = SentinelHubRequest(
        data_folder=data_folder,
        evalscript=evalscript,
        input_data=[
            SentinelHubRequest.input_data(
                data_collection=DataCollection.SENTINEL2_L1C,
                time_interval=time_interval,
                mosaicking_order='leastCC'
        )],
        responses=[
            SentinelHubRequest.output_response('default', MimeType.TIFF)
        ],
        bbox=loc_bbox,
        size=loc_size,
        config=config
    )
    
    outputs = request_all_bands.get_data(save_data=save_data)

    return outputs

I’d like to access the coordinates of each pixel of the outputs (or have two 1D arrays of latitudes and longitudes) so that I can point to the exact pixel I want. Is there a way to do so? It’s possible to create the matrix of coordinates using the Bbox but maybe there is a way to return it directly in the request. I could also of course save the output in a tif file, open it and get the coordinates but it’s very expensive when the operation is repeated many times.

Thank you in advance for your help!
Adrien

1 Like

Hi Adrien,

There are indeed some in-built functionalities to access coordinates in a sentinelhub-py request. For the sake of completeness, I will describe both ways (from and to pixel coordinates). My answer to your case is in point 2.

1. I would like to know the column/row number of a given coordinate

In this case, you can use the wgs84_to_pixel function found in geo_utils. This function takes a coordinate, and a transform vector to fetch the pixel coordinates of a given longitude/latitude pair (in WGS84, EPSG:4326). The little trick to using this function though is to specify the transform parameter in a UTM coordinate system, which you can fetch from your Bbox.

from sentinelhub import geo_utils

# Resolution in metres
resolution = 40
# Bounding box information and size of image
loc_bbox = BBox(bbox=coords, crs=CRS.WGS84)
loc_size = bbox_to_dimensions(loc_bbox, resolution=resolution)

# Convert bounding box to UTM coordinates
bb_utm = geo_utils.to_utm_bbox(loc_bbox)

# Fetch transform 
transf = bb_utm.get_transform_vector(resx=resolution, resy=resolution)

# And now you can query the pixel position of your WGS84 coordinates
pixel_pos = geo_utils.wgs84_to_pixel(125, 65, transf, utm_epsg=bb_utm.crs)

2. I would like to know the coordinates of a given row/column

In this case we can use the pixel_to_utm function from geo_utils. For a single coordinate pair, the approach is similar to point 1. But to get all the coordinates of the image, a little preparation is needed. We will make an array of pixel positions and then convert them to coordinates in UTM, before converting them to decimal degrees.

from sentinelhub import geo_utils

# Resolution in metres
resolution = 40
# Bounding box information and size of image
loc_bbox = BBox(bbox=coords, crs=CRS.WGS84)
loc_size = bbox_to_dimensions(loc_bbox, resolution=resolution)

# Convert bounding box to UTM coordinates
bb_utm = geo_utils.to_utm_bbox(loc_bbox)

# Fetch transform 
transf = bb_utm.get_transform_vector(resx=resolution, resy=resolution)

# Make an array of pixel values
# First the latitudes from the image size to a 2D grid
pix_lat = np.array(np.arange(0, loc_size[1]))
lats = np.array([pix_lat] * loc_size[0]).transpose()

# Second, the longitudes
pix_lon = np.array(np.arange(0, loc_size[0]))
lons = np.array([pix_lon] * loc_size[1])

# Convert the pixel positions to UTM
lon, lat = geo_utils.pixel_to_utm(lats, lons, transf)

# Convert your UTM pixel positions to WGS84 (EPSG:4326)
lon_degrees, lat_degrees = geo_utils.to_wgs84(lon, lat, bb_utm.crs)

Hope these code snippets help!

Maxim

Dear Maxim,

thank you a lot, it totally answers my question!
In order to reduce processing time, maybe I’d better directly extract the pixel at the targeted coordinates instead of downloading a larger area and then select it. Is it therefore possible to download data using a point instead of a bounding box? (in the same idea as your s3_extract tool). It doesnt seem possible to pass a point instead of a bounding box in a request. Again, I see a quick and dirty way to do it by giving the exact box size to get a 1x1 matrix as output but it is not the way you want to do it…!

thank you a lot,
Adrien

Hi Adrien,

You cannot specify a point in sentinelhub-py, but the approach that you suggest (a bounding box small enough to fall within 1 pixel) is a viable solution.

Although in my S3 Extract tool I got values using a point, you have to remember that the tool requires you to download the whole archive of images that you want to query.

This is where the power of Sentinel Hub truly shines: you will only need to download the values that you need.

As an example, imagine that you want to extract the values at 1 given pixel from a time series of 10 images (1 band): my tool requires you to download 10 images (10 * ~650 MB) then open the 10 bands = ~198978500 pixels; with sentinelhub-py you download 10 pixels worth of data.

Maxim

Hi there, I’ve been using this method and have been struggling to get it to work at high latitudes.

I suspect my problem is related to the lack of specification of the UTM zone in geo_utils.to_utm_bbox

But if anybody has any thoughts, please let me know