I have a list of large rasters that I would like to check systematically to see if they cover the same area as a list of points. Below is a program to walk the files and subdirectories of a folder and for each .tif, see if it contains at least 1 point. I used Google Colab to run this.
!pip install rasterio
!pip install geopandas
import rasterio
import rasterio.plot
import geopandas as gpd
from shapely.geometry import Point, Polygon
amuquall_df=gpd.read_file("/content/drive/MyDrive/amuquall/data/amuqall.shp")
print(amuquall_df.crs)
def plot_tif(tif_path, amuquall_df):
tif_file = rasterio.open(tif_path)
# print(tif_path)
# print(tif_file.shape)
# print(tif_file.bounds,tif_file.bounds[1])
# make a polygon of the tif
t_left = tif_file.bounds[0]
t_bottom = tif_file.bounds[1]
t_right = tif_file.bounds[2]
t_top = tif_file.bounds[3]
coords = [(t_left, t_bottom), (t_right, t_bottom), (t_right, t_top), (t_left, t_top)]
poly = Polygon(coords)
print("raster crs: ", tif_file.crs)
amuquall_df = amuquall_df.to_crs(tif_file.crs)
print("points crs: ", amuquall_df.crs)
tif_contains_points = False
for pt in amuquall_df['geometry']:
if(poly.contains(pt)):
tif_contains_points = True
break
return(tif_contains_points)
import os
usable_tifs = []
unusable_tifs = []
for dirpath, dirs, files in os.walk("/content/drive/MyDrive/crane_ane/Satellite_Data"):
for filename in files:
fname = os.path.join(dirpath,filename)
if fname.endswith('.tif'):
print(fname)
has_points = plot_tif(fname, amuquall_df)
if(has_points == True):
usable_tifs.append(fname)
if(has_points == False):
unusable_tifs.append(fname)
print("*** usable tifs ***")
print("\n".join(usable_tifs))
# print("*** UNusable tifs ***")
# print("\n".join(unusable_tifs))
Here is the output:
*** usable tifs *** /content/drive/MyDrive/crane_ane/Satellite_Data/SPOT/SPOT Images/Amuq.tif /content/drive/MyDrive/crane_ane/Satellite_Data/SPOT/SPOT Images/mosaic.tif /content/drive/MyDrive/crane_ane/Satellite_Data/SPOT/SPOT Images/nimadata1056733980.tif
… etc.