Geospatial Analysis: obtaining and pre-processing OpenSource satellite data

import os
import re
import time
from bs4 import BeautifulSoup
import numpy as np
import rasterio as rio
from osgeo import gdal
import requests
import dotenv
dotenv.load_dotenv(dotenv.find_dotenv())

# params fo LP DAAC
# you will require an EarthData login
host = fr'https://e4ftl01.cr.usgs.gov/MOTA/MCD64A1.061'
login = os.getenv('user')
password = os.getenv('pwd')

# list folders 
r = requests.get(host, verify=True, stream=True,auth=(login,password))
soup = BeautifulSoup(r.text, "html.parser")
folders = list()
for link in soup.findAll('a', attrs={'href': re.compile("\d{4}.\d{2}.\d{2}/")}):
    folders.append(link.get('href'))
print(f"{len(folders)} folders found")

# list files in folder
for f in folders[:1]:
    file_list = list()
    folder_url = f"{host}/{f}"
    r = requests.get(folder_url, verify=True, stream=True,auth=(login,password))
    soup = BeautifulSoup(r.text, "html.parser")
    for link in soup.findAll('a', attrs={'href': re.compile(".hdf$")}):
        file_list.append(link.get('href'))    
print(f"{len(file_list)} files found in folder n")

# USA tiles only
# use modis grid to slice only continents / areas of interest
# https://modis-land.gsfc.nasa.gov/MODLAND_grid.html
hreg = re.compile("h0[8-9]|h1[0-4]")
vreg = re.compile("v0[2-7]")
usa_files = list()
for fl in file_list:
    h = hreg.search(fl)
    if h:
        v = vreg.search(h.string)
        if v:
            usa_files.append(v.string)
print(f"{len(usa_files)} USA tiles found")

# download file from folder
# implement 25 file timeout to avoid rate limit

exceptions = list()
for f in folders[:1]:
    for e,fl in enumerate(usa_files[:]):
        if (e+1) % 10 == 0:
            print('sleep')
            time.sleep(5)
        else:
            print(fl,e+1)
            try:
                file_url = f"{host}/{f}/{fl}"
                r = requests.get(file_url, verify=True, stream=True,auth=(login,password))
                open(f'./data/raw/{fl}',"wb").write(r.content)
            except Exception as error:
                print(error)
                exceptions.append(fl)

for fl in usa_files[:]:       
    in_file = f"./data/raw/{fl}"
    out_file = f"./data/transformed/{fl.replace('.hdf','.tif')}"
    
    # open dataset
    dataset = gdal.Open(in_file,gdal.GA_ReadOnly)
    subdataset =  gdal.Open(dataset.GetSubDatasets()[0][0], gdal.GA_ReadOnly)
    
    # gdalwarp
    kwargs = {'format': 'GTiff', 'dstSRS': 'EPSG:4326'}
    ds = gdal.Warp(destNameOrDestDS=out_file,srcDSOrSrcDSTab=subdataset, **kwargs)
    del ds

# from rasterio.warp import calculate_default_transform, reproject, Resampling
# for fl in [fl for fl in usa_files if fl == 'MCD64A1.A2000306.h08v04.061.2021085165152.hdf']:       

#     file_name = f'./data/{fl}'

#     all_bands = []

#     with rio.open(file_name) as dataset:
#         # capture meta and CRS data
#         hdf4_meta = dataset.meta 
#         crs = dataset.read_crs()

#         # iterate data layers and select using name
#         for layer_name in dataset.subdatasets[:1]:
#             #print(layer_name)
#             with rio.open(layer_name) as subdataset:
#                 bounds = subdataset.bounds
#                 modis_meta = subdataset.profile
#                 all_bands.append(subdataset.read(1))

#         # prep metadata object
#         output_meta = modis_meta.copy()
#         output_meta['driver'] = 'GTiff'
#         output_meta['count'] = 1 # all_bands[0].shape[0]
        
#         with rio.open(out_path, "w", **kwargs) as dest:
#             dest.write(all_bands[0],indexes=1)
        
#         # reproject to 4326
#         dst_crs = "EPSG:4326"
#         transform, width, height = calculate_default_transform(output_meta.data['crs'],
#                                                                dst_crs,
#                                                                output_meta.data['width'],
#                                                                output_meta.data['height'],
#                                                                *bounds,
#                                                                dst_width=output_meta.data['width'],
#                                                                dst_height=output_meta.data['height'],)
#         kwargs = modis_meta.copy()
#         kwargs.update({'crs': dst_crs,
#                        'transform': transform,
#                        'width': width,
#                        'height': height})

#         # reproject and write to file e as a geotiff
#         out_path = f'{fl}.tif'
#         ino = all_bands[0]
#         ooo = all_bands[0].copy()
#         with rio.open(out_path, "w", **kwargs) as dest:
#             reproject(source=ino,
#                       destination=ooo,
#                       src_transform=modis_meta.data['transform'],
#                       src_crs=modis_meta.data['crs'],
#                       dst_transform=transform,
#                       dst_crs=dst_crs,
#                       resampling=Resampling.nearest)