1

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?

0
2

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('/tmp/test.tif') as src:

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

    # Generate a random window location that doesn't go outside the image
    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('/tmp/test_crop.tif', 'w', **profile) as dst:
        # Read the data from the window and write it to the output raster
        dst.write(src.read(window=window))
4
  • Worked beautifully. Thanks. – nrc Jul 16 '20 at 18:53
  • @user2856 that is great answer, I wonder if you couldplease elaborate how the generation of the random window works- not sure I understand the part of the src.width-xsize (and the same with the y) ? and why do we use random.randint? – Reut Nov 24 '20 at 9:38
  • @user2856 another thing i'm wondering is how do I do the connection between the coordinates to the pixel location ? – Reut Nov 24 '20 at 9:49
  • 1
    @Reut, the random.randint is because the OP wants random windows. The transform contains the geo coord info. – user2856 Nov 25 '20 at 20:10

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.