Using rasterio to crop image using pixel coordinates instead of geographic coordinates

Question

I'm trying to extract a series of random patches of image from a larger satellite image (from worldview3). I want to extract patches of uniform size (say 512x512 for example). If it were a normal image, any number of libraries could do this easily. But I want to use rasterio to retain the geographic information of the image patch.

So far the code I've written identifies the upper left pixel coordinate using

image.bounds

and using the resulting bounding box to get the height and the width. Then I use numpy to randomly generate a new random coordinate within the image extent:

new_col = np.random.randint(min(bounds.left, bounds.right), max(bounds.left, bounds.right)-(size+1))
new_row = np.random.randint(min(bounds.top, bounds.bottom), max(bounds.top, bounds.bottom)-(size+1))

And then I subtract the patch size (512) from the new row and column and use those new coordinates as the minx, miny, maxx, maxy and crop from there using

mask(image, shapes=coords, crop=True)

With a non-georeferenced image where the upper left is (0,0) and the lower right is (M,N) this works flawlessly. Similarly for a NAIP image this seems to work. But with the worldview3 image the size is not uniform. I'll get an image size like 700x2000 for example, where I want it to be 512x512. My thought was that the NAIP image I tested with was in UTM, but that the worldview3 image was in lat/long so that subtracting the patch size from the random row and column didn't translate to a uniform pixel size. I still think the error has to do with lat/long but after reprojecting the image to UTM the problem persists.

So, is there any way I can crop the worldview3 image using rasterio but using pixel coordinates instead of geographic coordinates so that I can crop a uniform image size from the larger image but still retain the geographic information of the cropped image patch?

Answer 1

You can use a window -rasterio.windows.Window to read by pixel offsets. The georeferencing can be easily calculated from the window using the source dataset window_transform method.

import random
import rasterio
from rasterio.windows import Window

with rasterio.open('input.tif') as src:

    # The size in pixels of your desired window
    xsize, ysize = 512, 512

    # Generate a random window origin (upper left) that ensures the window 
    # doesn't go outside the image. i.e. origin can only be between 
    # 0 and image width or height less the window width or height
    xmin, xmax = 0, src.width - xsize
    ymin, ymax = 0, src.height - ysize
    xoff, yoff = random.randint(xmin, xmax), random.randint(ymin, ymax)

    # Create a Window and calculate the transform from the source dataset    
    window = Window(xoff, yoff, xsize, ysize)
    transform = src.window_transform(window)

    # Create a new cropped raster to write to
    profile = src.profile
    profile.update({
        'height': xsize,
        'width': ysize,
        'transform': transform})

    with rasterio.open('output.tif', 'w', **profile) as dst:
        # Read the data from the window and write it to the output raster
        dst.write(src.read(window=window))