import pandas as pd
import geopandas as gpd
import folium
import os, shutil
from glob import glob

bounds = gpd.read_file('./data/processed/area_of_study_bounds.gpkg')

WRS_PATH = './data/external/Landsat8/wrs2_descending.zip'
LANDSAT_PATH = os.path.dirname(WRS_PATH)

!wget -P {LANDSAT_PATH} https://landsat.usgs.gov/sites/default/files/documents/wrs2_descending.zip

shutil.unpack_archive(WRS_PATH, os.path.join(LANDSAT_PATH, 'wrs2'))

wrs = gpd.GeoDataFrame.from_file('./data/external/Landsat8/wrs2/wrs2_descending.shp')

wrs.head()

wrs_intersection = wrs[wrs.intersects(bounds.geometry[0])]

paths, rows = wrs_intersection['PATH'].values, wrs_intersection['ROW'].values

Click on the Polygon to get the path and row value.

# Get the center of the map
xy = np.asarray(bounds.centroid[0].xy).squeeze()
center = list(xy[::-1])

# Select a zoom
zoom = 6

# Create the most basic OSM folium map
m = folium.Map(location=center, zoom_start=zoom, control_scale=True)

# Add the bounds GeoDataFrame in red
m.add_child(folium.GeoJson(bounds.__geo_interface__, name='Area of Study', 
                           style_function=lambda x: {'color': 'red', 'alpha': 0}))

# Iterate through each Polygon of paths and rows intersecting the area
for i, row in wrs_intersection.iterrows():
    # Create a string for the name containing the path and row of this Polygon
    name = 'path: %03d, row: %03d' % (row.PATH, row.ROW)
    # Create the folium geometry of this Polygon 
    g = folium.GeoJson(row.geometry.__geo_interface__, name=name)
    # Add a folium Popup object with the name string
    g.add_child(folium.Popup(name))
    # Add the object to the map
    g.add_to(m)

folium.LayerControl().add_to(m)
m.save('./images/10/wrs.html')
m

b = (paths > 23) & (paths < 26)
paths = paths[b]
rows = rows[b]

for i, (path, row) in enumerate(zip(paths, rows)):
    print('Image', i+1, ' - path:', path, 'row:', row)

Image 1  - path: 25 row: 27
Image 2  - path: 25 row: 28
Image 3  - path: 24 row: 28
Image 4  - path: 24 row: 27

s3_scenes = pd.read_csv('http://landsat-pds.s3.amazonaws.com/c1/L8/scene_list.gz', compression='gzip')
# google_scenes = pd.read_csv('https://storage.googleapis.com/gcp-public-data-landsat/index.csv.gz', compression='gzip')

s3_scenes.head(3)

# Empty list to add the images
bulk_list = []

# Iterate through paths and rows
for path, row in zip(paths, rows):

    print('Path:',path, 'Row:', row)

    # Filter the Landsat Amazon S3 table for images matching path, row, cloudcover and processing state.
    scenes = s3_scenes[(s3_scenes.path == path) & (s3_scenes.row == row) & 
                       (s3_scenes.cloudCover <= 5) & 
                       (~s3_scenes.productId.str.contains('_T2')) &
                       (~s3_scenes.productId.str.contains('_RT'))]
    print(' Found {} images\n'.format(len(scenes)))

    # If any scenes exists, select the one that have the minimum cloudCover.
    if len(scenes):
        scene = scenes.sort_values('cloudCover').iloc[0]

    # Add the selected scene to the bulk download list.
    bulk_list.append(scene)

Path: 25 Row: 27
 Found 2 images

Path: 25 Row: 28
 Found 6 images

Path: 24 Row: 28
 Found 3 images

Path: 24 Row: 27
 Found 6 images

bulk_frame = pd.concat(bulk_list, 1).T
bulk_frame

# Import requests and beautiful soup
import requests
from bs4 import BeautifulSoup

# For each row
for i, row in bulk_frame.iterrows():

    # Print some the product ID
    print('\n', 'EntityId:', row.productId, '\n')
    print(' Checking content: ', '\n')

    # Request the html text of the download_url from the amazon server. 
    # download_url example: https://landsat-pds.s3.amazonaws.com/c1/L8/139/045/LC08_L1TP_139045_20170304_20170316_01_T1/index.html
    response = requests.get(row.download_url)

    # If the response status code is fine (200)
    if response.status_code == 200:

        # Import the html to beautiful soup
        html = BeautifulSoup(response.content, 'html.parser')

        # Create the dir where we will put this image files.
        entity_dir = os.path.join(LANDSAT_PATH, row.productId)
        os.makedirs(entity_dir, exist_ok=True)

        # Second loop: for each band of this image that we find using the html <li> tag
        for li in html.find_all('li'):

            # Get the href tag
            file = li.find_next('a').get('href')

            print('  Downloading: {}'.format(file))

            # Download the files
            # code from: https://stackoverflow.com/a/18043472/5361345

            response = requests.get(row.download_url.replace('index.html', file), stream=True)

            with open(os.path.join(entity_dir, file), 'wb') as output:
                shutil.copyfileobj(response.raw, output)
            del response

EntityId: LC08_L1TP_025027_20170406_20170414_01_T1 

 Checking content:  

  Downloading: LC08_L1TP_025027_20170406_20170414_01_T1_B9.TIF.ovr
  Downloading: LC08_L1TP_025027_20170406_20170414_01_T1_B11.TIF
  Downloading: LC08_L1TP_025027_20170406_20170414_01_T1_B11.TIF.ovr
  ...
  ...
  ...

  Downloading: LC08_L1TP_025027_20170406_20170414_01_T1_B1.TIF.ovr

EntityId: LC08_L1TP_025028_20170913_20170928_01_T1 

  Checking content:  

  Downloading: LC08_L1TP_025028_20170913_20170928_01_T1_B4_wrk.IMD
  ...
  ...
  ...
  ...

import rasterio
from rasterio.transform import from_origin
from rasterio.warp import reproject, Resampling
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
%matplotlib inline

xmin, xmax, ymin, ymax = [], [], [], []

for image_path in glob(os.path.join(LANDSAT_PATH, '*/*B10.TIF')):
    with rasterio.open(image_path) as src_raster:
        xmin.append(src_raster.bounds.left)
        xmax.append(src_raster.bounds.right)        
        ymin.append(src_raster.bounds.bottom)        
        ymax.append(src_raster.bounds.top)        

fig, ax = plt.subplots(1, 1, figsize=(20, 15), subplot_kw={'projection': ccrs.UTM(16)})

ax.set_extent([min(xmin), max(xmax), min(ymin), max(ymax)], ccrs.UTM(16))

bounds.plot(ax=ax, transform=ccrs.PlateCarree())

for image_path in glob(os.path.join(LANDSAT_PATH, '*/*B10.TIF')):

    with rasterio.open(image_path) as src_raster:

        extent = [src_raster.bounds[i] for i in [0, 2, 1, 3]]  

        dst_transform = from_origin(src_raster.bounds.left, src_raster.bounds.top, 250, 250)

        width = np.ceil((src_raster.bounds.right - src_raster.bounds.left) / 250.).astype('uint')
        height = np.ceil((src_raster.bounds.top - src_raster.bounds.bottom) / 250.).astype('uint')

        dst = np.zeros((height, width))

        reproject(src_raster.read(1), dst, 
                  src_transform=src_raster.transform, 
                  dst_transform=dst_transform,
                  resampling=Resampling.nearest)

        ax.matshow(np.ma.masked_equal(dst, 0), extent=extent, transform=ccrs.UTM(16))


fig.savefig('./images/10/landsat.png')

shutil.rmtree('./data/external/Landsat8/LC08_L1TP_025027_20170406_20170414_01_T1')

scenes = s3_scenes[(s3_scenes.path == 25) & (s3_scenes.row == 27) & (s3_scenes.cloudCover <= 5)]
scenes

bulk_frame = pd.concat([bulk_frame.drop(679).T, scenes.loc[164949]], axis=1).T
bulk_frame