Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
NS-Nik opened this issue Apr 27, 2024 · 8 comments
Open

Comments

@NS-Nik
Copy link

NS-Nik commented Apr 27, 2024

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
Copy link
Contributor

mlubej commented May 7, 2024

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
Copy link
Author

NS-Nik commented May 21, 2024

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
Copy link
Contributor

mlubej commented May 21, 2024

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
Copy link
Author

NS-Nik commented May 21, 2024

@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
Copy link
Contributor

mlubej commented May 21, 2024

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
Copy link
Author

NS-Nik commented May 22, 2024

@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
Copy link
Contributor

mlubej commented May 27, 2024

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
Copy link
Author

NS-Nik commented Jun 10, 2024

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):
VCmELnLtl

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:

JpIdyPr2l

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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants