Convert pixel coordinates to latitude, longitude with the Sentinel-1 image

[NS-Nik]

Question
Hi!

I am trying to convert pixel coordinates to the geographical latitude, longitude. I am following steps, previously describe in this discussion: #51

However, I cannot retrieve information about CRS, coordinates of upper left corner and resolution.
I use rasterio:

src = rasterio.open(image)     
transform = src.transform 

But it seems like the standard GeoTIFF Sentinel-1 format does not contain the necessary information and transformation matrix looks like this:

print("Affine Transform:", transform)
Affine Transform:
 | 1.00, 0.00, 0.00|
| 0.00, 1.00, 0.00|
| 0.00, 0.00, 1.00|

Any ideas how to solve this problem?

Thank you a lot in advance for your help!

[mlubej]

Hi @NS-Nik!

Do you have gdalinfo in your command-line tool? Could you post the output for the tiff in question? Or share the file via some service, if it isn't too large.

[NS-Nik]

HI @mlubej!
Thank you for your reply!

The image size is quite big (~1.7 Gb), so here is the output from the gdalinfo:

metadata.txt

[mlubej]

Hi @NS-Nik!

Thanks for the command output. I wanted to check if the CRS info is actually available in the tif. Based on what you sent, it seems to be present, you just might not be accessing it correctly.

You mention that you are interested in the following:

However, I cannot retrieve information about CRS, coordinates of upper left corner and resolution.

I believe you should be able to extract the above by using other attributes rather than the transform, such as:

import rasterio

with rasterio.open(path_to_tif) as src:
    x1, y1, x2, y2 = src.bounds
    height = src.height
    width = src.width
    crs = src.crs.to_epsg()

res_x = (x2 - x1) // width
res_y = (y2 - y1) // height

This should work for the tiff based on what you sent me. You can then check if these values are also what you expect to get based on the gdalinfo output.

Hope this works!

[NS-Nik]

@mlubej Many thanks for your reply!

I still a little confused. Based on the information from the metadata (and from src.crs.to_epsg() command) , the srs is 4326 which is WGS 84.

However, the corner coordinates from the src.bounds do not look like they in the WGS 84.

Maybe you can take a look at the geotiff image?
https://drive.google.com/file/d/18MtufFSqfqm4F8xDG_FfJFwE3ewhShHh/view?usp=sharing

[mlubej]

Hmm, you are right @NS-Nik, the bounds or the CRS don't seem to match and there is some issue with the tif.

Is this the original tif, or did you get it somewhere else, or produced it yourself?

[NS-Nik]

@mlubej
This is the satellite image and it went through a series of processing steps after I downloaded it.
I have tried to do the same but with the different image (without implementing any additional processing steps).
Here is the metadata from this image:
meta.txt

The difference from the previous one - the information appears on the ground control points.

However, the problem stays the same: the Corner Coordinates are still not in the WGS 84 format.

[mlubej]

It seems that one or more of your processing steps after downloading are not done correctly, as the transformation information gets lost.

I would suggest you take a look into each of the processing steps after downloading and check when the data gets corrupted.

[NS-Nik]

Hi @mlubej!

Thank you for suggestion and sorry for the late response.

The problem is I need to batch process a hundreds of such images, which have already been downloaded, reprocessed and manually labeled. So, I don't have the opportunity to do the correct processing steps again and I am trying to find an approach to convert pixel coordinates to the lat, lon based on the information I have.

Here is how image looks when I just plot the data with imshow and plot one of the label with pixel coordinates (and this is how it should look):

Then I am trying to convert pixel coordinates and get the same result with lat, lon.

My code for the conversion:

#loading the data
data = gdal.Open(path).ReadAsArray()
latitude = data[2]
longitude = data[3]
hh = data[0]
hh[np.isnan(hh)] = 0

#the final png was flipped from right to left. Hence, I do the same with the initial data: 
hh = np.fliplr(hh)
latitude = np.fliplr(latitude)
longitude = np.fliplr(longitude)


#defining the latitude and longitude of the top left corner:
top_left_lat = latitude[0,0]
top_left_lon = longitude[0,0]
pixel_size_y = np.abs(latitude[0,0] - latitude[-1,0])/(latitude.shape[0]-1)
pixel_size_x = np.abs(longitude[0,-1]-longitude[0,0])/(latitude.shape[1]-1)


top_left_lat = latitude[0,0]
top_left_lon = longitude[0,0]
pixel_size_y = np.abs(latitude[0,0] - latitude[-1,0])/(latitude.shape[0]-1)
pixel_size_x = np.abs(longitude[0,-1]-longitude[0,0])/(latitude.shape[1]-1)


#convert yolo coordinates to pixel coordinates
x1 = 0.092308
y1 = 0.901066 
pixel_x  = x1 *  hh.shape[1]
pixel_y = y1 * hh.shape[0]

#define function for conversion

def pixel_to_geo(pixel_x, pixel_y):
     lon, lat = transform * (pixel_x, pixel_y)
     return lat, lon

#result
lat, lon = pixel_to_geo(pixel_x, pixel_y)

The final result looks like this:

As you can see, it misses a great deal. I also tried to calculate and add to the affine transformation the rotation angle, but result looks pretty much the same.

Any ideas how to solve this problem?

[YourVoicee]

Hi @NS-Nik!
Could you please clarify a few things in your code? I'm not sure what x1 and y1 stand for, and where their values are coming from. Also, in the "pixel_to_go" function, I don't understand what "transform" is. It's not being passed into the function, and it's not defined within the body of the function.