1

I'm trying to use rasterio or pyproj to convert latitude and longitudes into the SCREEN xy coordinates of a geotiff, so that I can superimpose 2D objects on it in a raster plotting library.

Specifically: I have a geotiff that is 512x512 in epsg:3857 representing a small area of satellite imagery.

I have a known lat long that is in that expanse somewhere. If it were in dead center, the transformation I want would yield 256,256 (50%,50% of 512,512).

How to do it? I tried:

import rasterio as rio
latitude, longitude = 39.60378691188893, -104.9429961240613
dataset = rio.open('satellite_downloads/geotiff/home.tif', crs='epsg:3857')
pprint(dataset.crs)
pprint(dataset.transform)
py, px = dataset.index(x= longitude, y=latitude)
print('Pixel Y, X coords: {}, {}'.format(py, px))

But this returns

CRS({'init': 'epsg:3857'}) Affine(0.14929107090938487, 0.0, -11682239.10936522, 0.0, -0.1492910709366697, 4808569.838743366)

Pixel Y, X coords: 32209094, 782507220

How do I convert that into a value between 0 - 511 (my raster xy range)?

4
  • Have you tried swapping lat/lon to index(x=lon, y=lat)? – mikewatt Apr 29 at 17:20
  • hi! yes, thank you, that was my mistake, but doesn't reproject the coords (updated ticket with that fix). Any ideas how do I convert 32209094, 782507220 into my screen pixel dimensions? – metalaureate Apr 29 at 17:27
  • Are your coordinates also in web mercator? If they're in WGS84 or something else then you'll need to reproject them to match the raster (or the other way around) – mikewatt Apr 29 at 19:19
  • The coordinates are latitude and longitude. If I convert them Mercator, I get y=-11682195.3, x=4808529.7, which produces 110460491, 110460516. None of that can be right. Conceptually, I am missing the transformation from EPSG:3857 to the pixels in my raster image (512x512). I don't know how to do that transformation. – metalaureate Apr 29 at 19:57
1

I was overthinking this, which may be why no one has posted an answer before:

We just need to calculate the x,y proportion of the image bounding box and multiply the pixel size of the satellite image to get the x,y of the image.

from pprint import pprint
from pyproj import Proj, transform
import rasterio as rio

def image_latlon_pxpy(latitude, longitude):
    dataset = rio.open('satellite_downloads/geotiff/home.tif', crs='epsg:3857')
    coords = transform(Proj(init='epsg:4326'), Proj(init='epsg:3857'), longitude, latitude)
    px, py = coords[0], coords[1]
    print(px,py)
    px_pc = (px - dataset.bounds.left) / (dataset.bounds.right - dataset.bounds.left)
    py_pc = (dataset.bounds.top - py) / (dataset.bounds.top - dataset.bounds.bottom)
    return (px_pc*dataset.width, py_pc*dataset.height)

pprint(image_latlon_pxpy(39.603786, -104.9429))
(256,256)

I assume this may be increasingly inaccurate for large satellite image tiles without spherical projection, but for small tiles it works fine.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.