Checking a list of rasters to see if they correspond to a list of points.

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


def plot_tif(tif_path, amuquall_df):
    tif_file =
    # 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: ",
    amuquall_df = amuquall_df.to_crs(
    print("points crs: ",
    tif_contains_points = False
    for pt in amuquall_df['geometry']:
            tif_contains_points = True

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'):
            has_points = plot_tif(fname, amuquall_df)
            if(has_points == True):
            if(has_points == False):
print("*** 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.

Plotting Shape files using Google Colab and Earth Engine

I like to use QGIS to examine shapefiles, but sometimes I use a computer that does not have it installed. Below are some notes about plotting shapefiles using web interfaces I can access on most computers: Google Colab and Earth Engine.

Ashwani Dhankhar’s article “Mapping with Matplotlib, Pandas, Geopandas and Basemap in Python” shows how to plot a shapefile on a basemap in Google Colab. Below is code in Google Colab to do this. Following that, is a note about Google Earth Engine.

# !apt-get install libgeos-3.6.2 (I am not sure if this is needed)
!apt-get install libgeos-dev
!pip install
!pip install geopandas
!pip install contextily

Restart the runtime.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import geopandas as gpd
import pandas as pd
import contextily as ctx

The shapefile I want to view (amuqall.shp) is on my Google drive. I wanted to view its geographic extent on a basemap.

# Change crs to one compatible with basemap
amuquall_df = amuquall_df.to_crs(epsg=3857)
ax = amuquall_df.plot(figsize=(20, 20), alpha=0.5, edgecolor='k')