Hello everyone,

In the picture below I am using rasterio to open a .tif file and then transform the those bounds to get the upper left and lower right points of the dataset. My question is around how to convert a point that is not from the bounds and is just an x,y point in meters in a specific epsg. I want to be able to convert multiple points that are in the x,y meter format to a lat/long format all within this epsg.

Another question I have is around how applicable this is to creating new distances. For example is I am at the upper left spot and increase the x by 5 meters and convert that coordinate back to lat/long, does that mean that the two lat/long coordinates are now exactly 5 meters apart? 

I am quite new to the GIS world so any help would be greatly appreciated, thanks!

I 

Hi,

From the screenshot you have provided, we can see that the CRS of your image (EPSG:32617) is a UTM zone, so its unit is indeed in meters. That means that if you take a point with coordinates expressed in the CRS of your image, that you add 5 meters to the x coordinates, and that you convert both points to lat/lon, your two points will indeed be separated by 5 meters.

Here is a code snippet showing what you want to achieve, that uses the pyproj package to do the coordinate conversion:

import rasterio
import pyproj

with rasterio.open("S2A_MSIL1C_20190701T174911_N0207_R141_T13UFQ_20190701T212910__B11.jp2") as src:
    bounds = src.bounds
    crs = src.crs

left, bottom, right, top = bounds

in_proj = pyproj.Proj(init="epsg:{}".format(crs.to_epsg()))
out_proj = pyproj.Proj(init="epsg:4326")

lons, lats = pyproj.transform(in_proj, out_proj, [left, left + 5], [top, top])

for lat, lon in zip(lats, lons):
    print(lat, lon)

This snippet uses a Sentinel 2 band with a different CRS than yours, but it should work the same way.

A few things to note:

• you can directly access the left, bottom, right and top bounds from the bounds object, there is no need to use the transform for that.
• EPSG:4326 is the CRS for coordinates expressed in lat/lon
• pyproj.transform can take either single x, y coordinates, or iterables of x, y coordinates (as in the example above where I pass a list of two x's and two y's). If you want to convert a lot of points between two CRS's, it is much more efficient to pass them as a list in a single call to pyproj.transform because the transformation is vectorized.

Guillaume Lostis

Guillaume, thanks for the swift reply, this sounds exactly like what I'm looking for. I will work on implementing it into my code now, thank again :)