How to get mask from a huge .tif file in rasterio python?

user_4921122 0
Viewed 1302

I have a large (~3GB) .tif file with Geo information embedded in it.
I have a few Polygons (which are expressed by GPS coordinates) saved in a .shp file that describe specific areas in the image that interest me.

I want to get a crop of the image where each polygon refers to, and also the mask of the polygon on that crop.

I can read a window based on the bounding box from the polygons,
but then I can't match the polygon inside the numpy array created from the rasterio.read function:

from math import ceil
import rasterio
import fiona

with fiona.open("shapes.shp", "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]

tif_fn = 'large_file.tif'
my_tif = rasterio.open(tif_fn)

bound_1 = rasterio.features.bounds(shapes[0])
bbox_1 = rasterio.windows.from_bounds(*bound_1, my_tif.transform)
window_transform1 = rasterio.windows.transform(window=bbox_1, transform=my_tif.transform)
mask = rasterio.features.geometry_mask([shapes[0]],out_shape=(ceil(bbox_1.width), ceil(bbox_1.height)), 
                                       transform=window_transform, invert=True)
img_crop = my_tif.read([1,2,3], window=bbox_1) # pretty fast, ~2 seconds

plt.imshow(img_crop)
plt.imshow(mask,alpha=0.2)
plt.show() # bad match of image and mask...

I've tried using rasterio.mask.mask with crop=True based on this tutorial, but it takes way too long on huge files. (~50 seconds)

# takes 50 seconds...
out_image, out_transform = rasterio.mask.mask(my_tif, shapes, crop=True, filled=True) 
out_meta = rsrc.meta

Is there a way to maybe to make a sub-Dataset object and get a mask from that?
Or a way to get the crop and mask together?

thanks!

rasterio python
1 Answers

first option is simply fixed, just use the cropped image shape for the mask:

img_crop = my_tif.read([1,2,3], window=bbox_1) # pretty fast, ~2 seconds
mask = rasterio.features.geometry_mask([shapes[0]],out_shape=(img_crop.shape[1], img_crop.shape[2]), 
                                       transform=window_transform1, invert=True)

the problem with the rasterio.mask.mask was that I was trying to make a mask for all my shapes, instead of doing it for each one. I needed to change the code like this:

# takes 6 seconds, way better.
out_image, out_transform = rasterio.mask.mask(my_tif, [shapes[0]], crop=True, filled=True) 
out_meta = rsrc.meta
original_image_cropped = my_tif.read([1,2,3], window=bbox_1)
plt.imshow(original_image[[0,1,2]])
plt.imshow(out_image[0], alpha=0.2)
plt.show()

# out_image is the crop, and read

notice I only mask one shape, not all of them.

user_49211220
Related