I'm not a GIS person. I haven't found or at least identified a solution to my problem, so decided to ask.

I have a big rasters and a list of thousands of coordinates. I need to find the pixel location corresponding to a given lat/lon pair in the raster, and then clip an NxN rectangle with the pixel located in the middle of it and save the clip on disk. I haven't for the life of me figured out yet how I could find the pixel location using coordinates, and I've been on it a few days now. I found that with snappy I should be able to do something like that, but unfortunately it doesn't run on my system (anaconda3 with Python 3.6).

How can I go about solving this problem?

  • Going from coordinates in a given CRS to array indices and vice-versa is what the affine library does. Any rasterio dataset has an Affine which can be accessed via the transform attribute – Loïc Dutrieux Oct 22 '18 at 21:50
  • You have a tag for GDAL, is that the API you're interested in? In GDAL you open your raster and retrieve the GeoTransform which is 6 double values see the tutorial gdal.org/gdal_tutorial.html, column is ( X - GeoTransform[0] ) / GeoTransform[1] and row is ( GeoTransform[3] - Y ) / GeoTransform[5] because the reference is top left and Y cell sizes are negative. – Michael Stimson Oct 22 '18 at 22:45
  • That formula only works if the raster is georeferenced docs.qgis.org/2.18/en/docs/training_manual/forestry/… and is in the same spatial reference as your input points. – Michael Stimson Oct 22 '18 at 22:57
  • Sorry, should've been more explicit about the APIs. I'm interested in both GDAL and rasterio, although I prefer the latter. In any case, Luke's answer is exactly what I've been trying to do, and it's with rasterio. Thanks anyways for your help! – GISjoe Oct 23 '18 at 16:33

A rasterio way of doing this is pretty simple. Note this requires your raster be in the same projection as your coordinates. You can of course project your coordinates on the fly, but that's another question...

import rasterio as rio

infile = r"C:\Temp\test.tif"
outfile = r'C:\Temp\test_{}.tif'
coordinates = (
    (130.5, -25.5) , # lon, lat of ~centre of Australia
    (146.0, -42.0) , # lon, lat of ~centre of Tasmania

# Your NxN window
N = 3

# Open the raster
with rio.open(infile) as dataset:

    # Loop through your list of coords
    for i, (lon, lat) in enumerate(coordinates):

        # Get pixel coordinates from map coordinates
        py, px = dataset.index(lon, lat)
        print('Pixel Y, X coords: {}, {}'.format(py, px))

        # Build an NxN window
        window = rio.windows.Window(px - N//2, py - N//2, N, N)

        # Read the data in the window
        # clip is a nbands * N * N numpy array
        clip = dataset.read(window=window)

        # You can then write out a new file
        meta = dataset.meta
        meta['width'], meta['height'] = N, N
        meta['transform'] = rio.windows.transform(window, dataset.transform)

        with rio.open(outfile.format(i), 'w', **meta) as dst:
  • Thanks, Luke, your answer is exactly what I've been trying to do with no success. I just tried it with the raster and coordinates I have, and it seems to be working like a charm. – GISjoe Oct 23 '18 at 16:35

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.