Determine Valid MGRS Sqaures to Process¶

In [1]:

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
import pystac
import pystac_client
import stackstac
import matplotlib.pyplot as plt
import xarray as xr
from datetime import datetime
import matplotlib
import ipyleaflet
import sys
import os
import dask_gateway
import planetary_computer
from rechunker import rechunk
sys.path.append('../sar_snowmelt_timing')
import s1_rtc_bs_utils
import contextily as ctx
import rioxarray as rxr
import pathlib
import glob
import re
import time
import fsspec

In [2]:

cluster = dask_gateway.GatewayCluster(shutdown_on_close=False)
client = cluster.get_client()
cluster.adapt(minimum=90, maximum=100)
print(client.dashboard_link)

https://pccompute.westeurope.cloudapp.azure.com/compute/services/dask-gateway/clusters/prod.d3da69dfe3f140358199cdd825157e7e/status

In [3]:

gj = '../input/shapefiles/western_us.geojson'
bbox_gdf = gpd.read_file(gj)

In [4]:

url = 'https://github.com/scottyhq/mgrs/raw/main/MGRS_LAND.parquet' # Scott created an MGRS parquet file here https://github.com/scottyhq/mgrs
with fsspec.open(url) as file:
    mgrs_gdf = gpd.read_parquet(file)

In [5]:

mgrs_gdf_clipped = mgrs_gdf.clip(bbox_gdf)

In [6]:

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
    )

bbox = mgrs_gdf_clipped.total_bounds

search = catalog.search(
    collections=["esa-worldcover"],
    bbox=bbox,
    )

items = list(search.get_items())
items

stack_lc = stackstac.stack(items, bounds_latlon=bbox,resolution=0.00833333)#ts_ds.resolution or ts_ds.rio.resolution()[0]
worldcover = stack_lc.min(dim='time').squeeze()

In [7]:

snow_classes = [ # https://zenodo.org/record/2626737#.ZDMMf3bMIQ8
#    0, #Little-to-no snow
    1, #Indeterminate due to clouds
    2, #Ephemeral snow
    3, #Seasonal snow
]

classes = [ # page 13 of https://esa-worldcover.s3.amazonaws.com/v100/2020/docs/WorldCover_PUM_V1.0.pdf
#    10, # treecover
    20, # shrubland
    30, # grassland
    40, # cropland
#    50, # built-up
    60, #bare / sparse vegetation
    70, # snow and ice
#    80, # permanent water bodies
    90, # herbaceous wetlands
    95, # mangroves
    100 # loss and lichen
]

In [8]:

worldcover = worldcover.rio.write_crs(rio.CRS.from_epsg(4326))

In [ ]:

snow_mask = rxr.open_rasterio('../input/SnowClass/global_MODIS_snow_classes_byte.tif')
    
snow_mask_clip = snow_mask.sel(x=slice(bbox[0],bbox[2]),y=slice(bbox[3],bbox[1])).squeeze()  
snow_mask_proj = snow_mask_clip.rio.reproject_match(worldcover)

In [ ]:

valid_pixels = snow_mask.isin(snow_classes) & worldcover.isin(classes)

In [12]:

f,ax=plt.subplots(figsize=(10,10))
valid_pixels.plot(ax=ax,vmax=1)
mgrs_gdf_clipped.boundary.plot(ax=ax,color='black')
ctx.add_basemap(ax=ax, crs=valid_pixels.rio.crs, source=ctx.providers.Stamen.Terrain,attribution=False)

In [14]:

pixels = []
for i,row in mgrs_gdf_clipped.iterrows():
    num_pixels = snow_present.rio.clip([row.geometry])>0
    pixels.append(num_pixels.where(num_pixels>0).count())

In [15]:

mgrs_gdf_clipped['bare_ground_and_snow_pixels'] = np.array(pixels)

In [16]:

f,ax=plt.subplots(figsize=(10,10))
mgrs_gdf_clipped.plot(ax=ax,column='bare_ground_and_snow_pixels',)
ctx.add_basemap(ax=ax, crs=valid_pixels.rio.crs, source=ctx.providers.Stamen.Terrain,attribution=False)

In [17]:

tiles = mgrs_gdf_clipped.loc[mgrs_gdf_clipped['bare_ground_and_snow_pixels']>50].sort_values(by='bare_ground_and_snow_pixels',ascending=False).tile

In [25]:

file = open('../input/MGRS_square_list/westernUS_square_list.txt','w')
for item in tiles:
    file.write(item+"\n")
file.close()

In [ ]: