In [1]:

#!mamba env create -f '../../environment.yml'
#!mamba env update -f '../../environment.yml' --prune

#!mamba env create -f '../../environment.yml' #!mamba env update -f '../../environment.yml' --prune

Testing examples from remote_sensing¶

In [31]:

# this block is for developing the module, comment out when using the module, and uncomment import easysnowdata
%load_ext autoreload
%autoreload 2
%aimport easysnowdata

# this block is for developing the module, comment out when using the module, and uncomment import easysnowdata %load_ext autoreload %autoreload 2 %aimport easysnowdata

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

In [32]:

#import easysnowdata
import geopandas as gpd
import rioxarray as rxr
import xarray as xr
import shapely
import dask
import pystac_client
import planetary_computer
import odc.stac
import matplotlib.pyplot as plt
import matplotlib.colors
import datetime
import pandas as pd
import numpy as np
import contextily as ctx
import rasterio as rio

import datetime

today = datetime.datetime.now().strftime("%Y-%m-%d")

#import easysnowdata import geopandas as gpd import rioxarray as rxr import xarray as xr import shapely import dask import pystac_client import planetary_computer import odc.stac import matplotlib.pyplot as plt import matplotlib.colors import datetime import pandas as pd import numpy as np import contextily as ctx import rasterio as rio import datetime today = datetime.datetime.now().strftime("%Y-%m-%d")

In [33]:

bbox_gdf = gpd.read_file(
    "https://github.com/egagli/easysnowdata/raw/main/docs/examples/mt_rainier.geojson"
)

bbox_gdf = gpd.read_file( "https://github.com/egagli/easysnowdata/raw/main/docs/examples/mt_rainier.geojson" )

Copernicus Global Land Service forest cover fraction¶

In [34]:

forest_cover_fraction_da = easysnowdata.remote_sensing.get_forest_cover_fraction(bbox_gdf)
#forest_cover_fraction_da

forest_cover_fraction_da = easysnowdata.remote_sensing.get_forest_cover_fraction(bbox_gdf) #forest_cover_fraction_da

In [35]:

f,ax = forest_cover_fraction_da.attrs['example_plot'](forest_cover_fraction_da)

f,ax = forest_cover_fraction_da.attrs['example_plot'](forest_cover_fraction_da)

Sturm & Liston 2021 seasonal snow classification¶

In [6]:

snow_classification_da = easysnowdata.remote_sensing.get_seasonal_snow_classification(bbox_gdf)
#snow_classification_da

snow_classification_da = easysnowdata.remote_sensing.get_seasonal_snow_classification(bbox_gdf) #snow_classification_da

In [7]:

f,ax = snow_classification_da.attrs['example_plot'](snow_classification_da)

f,ax = snow_classification_da.attrs['example_plot'](snow_classification_da)

¶

Wrzesien et al 2019 global seasonal mountain snow mask¶

In [36]:

mountain_snow_da = easysnowdata.remote_sensing.get_seasonal_mountain_snow_mask(bbox_gdf)
#mountain_snow_da

mountain_snow_da = easysnowdata.remote_sensing.get_seasonal_mountain_snow_mask(bbox_gdf) #mountain_snow_da

This function takes a moment, getting zipped file from zenodo...

In [37]:

f, ax = mountain_snow_da.attrs['example_plot'](mountain_snow_da)

f, ax = mountain_snow_da.attrs['example_plot'](mountain_snow_da)

ESA WorldCover¶

In [10]:

worldcover_da = easysnowdata.remote_sensing.get_esa_worldcover(bbox_gdf)
#worldcover_da

worldcover_da = easysnowdata.remote_sensing.get_esa_worldcover(bbox_gdf) #worldcover_da

In [11]:

f, ax = worldcover_da.attrs['example_plot'](worldcover_da)

f, ax = worldcover_da.attrs['example_plot'](worldcover_da)

Sentinel-2¶

In [12]:

s2 = easysnowdata.remote_sensing.Sentinel2(
    bbox_input=bbox_gdf,
    start_date="2024-05-31",
    end_date="2024-07-07",
    catalog_choice="planetarycomputer", # can also choose "earthsearch". one might be faster than the other depending on where you are running the code
    resolution=80,
)
#s2.data

s2 = easysnowdata.remote_sensing.Sentinel2( bbox_input=bbox_gdf, start_date="2024-05-31", end_date="2024-07-07", catalog_choice="planetarycomputer", # can also choose "earthsearch". one might be faster than the other depending on where you are running the code resolution=80, ) #s2.data

Data searched. Access the returned seach with the .search attribute.
Data retrieved. Access with the .data attribute. Data CRS: WGS 84 / UTM zone 10N.
Nodata values removed from the data. In doing so, all bands converted to float32. To turn this behavior off, set remove_nodata=False.
Data acquired after January 25th, 2022 harmonized to old baseline. To override this behavior, set harmonize_to_old=False.
Data scaled to float reflectance. To turn this behavior off, set scale_data=False.
Metadata retrieved. Access with the .metadata attribute.

Get Sentinel-2 RGB imagery¶

In [13]:

s2.get_rgb()
#s2.rgb
#s2.rgba
#s2.rgb_percentile
#s2.rgb_clahe

s2.get_rgb() #s2.rgb #s2.rgba #s2.rgb_percentile #s2.rgb_clahe

/home/eric/miniconda3/envs/easysnowdata/lib/python3.10/site-packages/odc/geo/_rgba.py:56: RuntimeWarning: invalid value encountered in cast
  return x.astype("uint8")

/home/eric/miniconda3/envs/easysnowdata/lib/python3.10/site-packages/odc/geo/_rgba.py:56: RuntimeWarning: invalid value encountered in cast
  return x.astype("uint8")

RGB data retrieved.
Access with the following attributes:
.rgb for raw RGB,
.rgba for RGBA,
.rgb_percentile for percentile RGB,
.rgb_clahe for CLAHE RGB.
You can pass in percentile_kwargs and clahe_kwargs to adjust RGB calculations, check documentation for options.

In [14]:

f = s2.rgba.plot.imshow(col='time',col_wrap=5, robust=False)

for ax, time, in zip(f.axs.flat, s2.rgba.time.values):
    local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles')
    ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}')
    ax.axis('off')
    ax.set_aspect('equal')
    
f.fig.tight_layout()
f.fig.dpi = 300
f.fig.suptitle('Sentinel-2 RGBA', fontsize=16, y=1.02)

f = s2.rgba.plot.imshow(col='time',col_wrap=5, robust=False) for ax, time, in zip(f.axs.flat, s2.rgba.time.values): local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles') ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}') ax.axis('off') ax.set_aspect('equal') f.fig.tight_layout() f.fig.dpi = 300 f.fig.suptitle('Sentinel-2 RGBA', fontsize=16, y=1.02)

/home/eric/miniconda3/envs/easysnowdata/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(

Out[14]:

Text(0.5, 1.02, 'Sentinel-2 RGBA')

In [15]:

f = s2.rgb_clahe.plot.imshow(col='time',col_wrap=5, robust=False)

for ax, time, in zip(f.axs.flat, s2.rgb.time.values):
    local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles')
    ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}')
    ax.axis('off')
    ax.set_aspect('equal')
    
f.fig.tight_layout()
f.fig.dpi = 300
f.fig.suptitle('Sentinel-2 RGB w/ clahe equalization', fontsize=16, y=1.02)

f = s2.rgb_clahe.plot.imshow(col='time',col_wrap=5, robust=False) for ax, time, in zip(f.axs.flat, s2.rgb.time.values): local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles') ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}') ax.axis('off') ax.set_aspect('equal') f.fig.tight_layout() f.fig.dpi = 300 f.fig.suptitle('Sentinel-2 RGB w/ clahe equalization', fontsize=16, y=1.02)

Out[15]:

Text(0.5, 1.02, 'Sentinel-2 RGB w/ clahe equalization')

In [16]:

f,axs = plt.subplots(1,3,figsize=(12,7))

s2.rgba.isel(time=0).plot.imshow(ax=axs[0])
axs[0].set_title('RGBA')
s2.rgb_percentile.isel(time=0).plot.imshow(ax=axs[1])
axs[1].set_title('RGB w/ percentile stretch')
s2.rgb_clahe.isel(time=0).plot.imshow(ax=axs[2])
axs[2].set_title('RGB w/ clahe equalization')

for ax in axs:
    ax.axis('off')
    ax.set_aspect('equal')

f.tight_layout()

f,axs = plt.subplots(1,3,figsize=(12,7)) s2.rgba.isel(time=0).plot.imshow(ax=axs[0]) axs[0].set_title('RGBA') s2.rgb_percentile.isel(time=0).plot.imshow(ax=axs[1]) axs[1].set_title('RGB w/ percentile stretch') s2.rgb_clahe.isel(time=0).plot.imshow(ax=axs[2]) axs[2].set_title('RGB w/ clahe equalization') for ax in axs: ax.axis('off') ax.set_aspect('equal') f.tight_layout()

Get Sentinel-2 band indicies (NDSI, NDVI, NDWI, EVI, NDBI implemented)¶

In [17]:

s2.get_ndvi()
#s2.ndvi

s2.get_ndvi() #s2.ndvi

NDVI data calculated. Access with the .ndvi attribute.

In [18]:

f = s2.ndvi.plot.imshow(col='time',col_wrap=5)

for ax, time, in zip(f.axs.flat, s2.ndvi.time.values):
    local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles')
    ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}')
    ax.axis('off')
    ax.set_aspect('equal')

f.fig.suptitle('Sentinel-2 NDVI time series',fontsize=16,y=1.02)

f = s2.ndvi.plot.imshow(col='time',col_wrap=5) for ax, time, in zip(f.axs.flat, s2.ndvi.time.values): local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles') ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}') ax.axis('off') ax.set_aspect('equal') f.fig.suptitle('Sentinel-2 NDVI time series',fontsize=16,y=1.02)

Out[18]:

Text(0.5, 1.02, 'Sentinel-2 NDVI time series')

Clouds get in the way! Let's check out the scene classification layer (SCL)¶

In [19]:

#s2.data.scl

#s2.data.scl

In [20]:

f, ax = s2.data.scl.attrs['example_plot'](s2.data.scl.sel(time='2024-05-31',method='nearest'))

f, ax = s2.data.scl.attrs['example_plot'](s2.data.scl.sel(time='2024-05-31',method='nearest'))

In [21]:

f, ax = s2.data.scl.attrs['example_plot'](s2.data.scl)

f, ax = s2.data.scl.attrs['example_plot'](s2.data.scl)

Let's say we want to view cloud free images or calculate NDSI while mitigating clouds. We can do that by masking our data based on the SCL. Let's mask the data we are not interested in. You can mask data by passing boolean values to .mask_data(). By default, no data, saturated pixels, topographic shadows, cloud shadows, cloud medium probability, cloud high probability, and thin cirrus are all removed.¶

In [22]:

s2.mask_data()

s2.mask_data()

Removed pixels with the following scene classification values:
No Data (Missing data)
Saturated or defective pixel
Topographic casted shadows
Cloud shadows
Cloud medium probability
Cloud high probability
Thin cirrus

In [23]:

# get rgb again, now with masked data
s2.get_rgb()
#s2.rgb_clahe

# get rgb again, now with masked data s2.get_rgb() #s2.rgb_clahe

RGB data retrieved.
Access with the following attributes:
.rgb for raw RGB,
.rgba for RGBA,
.rgb_percentile for percentile RGB,
.rgb_clahe for CLAHE RGB.
You can pass in percentile_kwargs and clahe_kwargs to adjust RGB calculations, check documentation for options.

In [24]:

f = s2.rgb_clahe.plot.imshow(col='time',col_wrap=5, robust=True)

for ax, time, in zip(f.axs.flat, s2.rgb.time.values):
    local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles')
    ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}')
    ax.axis('off')
    ax.set_aspect('equal')
    
f.fig.tight_layout()
f.fig.dpi = 300
f.fig.suptitle('Sentinel-2 RGB w/ clahe equalization', fontsize=16, y=1.02)

f = s2.rgb_clahe.plot.imshow(col='time',col_wrap=5, robust=True) for ax, time, in zip(f.axs.flat, s2.rgb.time.values): local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles') ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}') ax.axis('off') ax.set_aspect('equal') f.fig.tight_layout() f.fig.dpi = 300 f.fig.suptitle('Sentinel-2 RGB w/ clahe equalization', fontsize=16, y=1.02)

Out[24]:

Text(0.5, 1.02, 'Sentinel-2 RGB w/ clahe equalization')

In [25]:

s2.get_ndsi()
#s2.ndsi

s2.get_ndsi() #s2.ndsi

NDSI data calculated. Access with the .ndsi attribute.

In [26]:

f = s2.ndsi.plot.imshow(col='time',col_wrap=5)

for ax, time, in zip(f.axs.flat, s2.ndsi.time.values):
    local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles')
    ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}')
    ax.axis('off')
    ax.set_aspect('equal')

f.fig.dpi = 300
f.fig.suptitle('Sentinel-2 NDSI time series',fontsize=16,y=1.02)

f = s2.ndsi.plot.imshow(col='time',col_wrap=5) for ax, time, in zip(f.axs.flat, s2.ndsi.time.values): local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles') ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}') ax.axis('off') ax.set_aspect('equal') f.fig.dpi = 300 f.fig.suptitle('Sentinel-2 NDSI time series',fontsize=16,y=1.02)

Out[26]:

Text(0.5, 1.02, 'Sentinel-2 NDSI time series')

Sentinel-1¶

In [27]:

s1 = easysnowdata.remote_sensing.Sentinel1(
    bbox_input=bbox_gdf, start_date="2022-07-01", end_date="2022-07-31", resolution=80
)
#s1.data

s1 = easysnowdata.remote_sensing.Sentinel1( bbox_input=bbox_gdf, start_date="2022-07-01", end_date="2022-07-31", resolution=80 ) #s1.data

Data searched. Access the returned seach with the .search attribute.
Data retrieved. Access with the .data attribute. Data CRS: WGS 84 / UTM zone 10N.
Metadata retrieved. Access with the .metadata attribute.
Falsely low scenes and border noise removed from the data.
Added relative orbit number and orbit state as coordinates to the data.
Linear power units converted to dB. Convert back to linear power units using the .db_to_linear() method.

In [28]:

f = s1.data['vv'].plot.imshow(col='time',col_wrap=5, vmin=-15, vmax=2, cmap='gray')

for ax, time, in zip(f.axs.flat, s1.data['vv'].time.values):
    local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles')
    ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}')
    ax.axis('off')
    ax.set_aspect('equal')

f.fig.suptitle('Sentinel-1 RTC backscatter time series',fontsize=16,y=1.02)

f = s1.data['vv'].plot.imshow(col='time',col_wrap=5, vmin=-15, vmax=2, cmap='gray') for ax, time, in zip(f.axs.flat, s1.data['vv'].time.values): local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles') ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}') ax.axis('off') ax.set_aspect('equal') f.fig.suptitle('Sentinel-1 RTC backscatter time series',fontsize=16,y=1.02)

Out[28]:

Text(0.5, 1.02, 'Sentinel-1 RTC backscatter time series')

In [29]:

f, ax = plt.subplots(figsize=(10, 10))

s1.metadata.plot(
    "sat:relative_orbit",
    ax=ax,
    edgecolor="black",
    categorical=True,
    legend=True,
    alpha=0.2,
)
bbox_gdf.plot(ax=ax, edgecolor="black", facecolor="none")
ctx.add_basemap(ax, crs=s1.metadata.crs, source=ctx.providers.Esri.WorldImagery)

f, ax = plt.subplots(figsize=(10, 10)) s1.metadata.plot( "sat:relative_orbit", ax=ax, edgecolor="black", categorical=True, legend=True, alpha=0.2, ) bbox_gdf.plot(ax=ax, edgecolor="black", facecolor="none") ctx.add_basemap(ax, crs=s1.metadata.crs, source=ctx.providers.Esri.WorldImagery)

Harmonized Landsat Sentinel-2 v2.0¶

In [30]:

# !CPL_VSIL_CURL_USE_HEAD=FALSE
# !GDAL_DISABLE_READDIR_ON_OPEN=YES
# !GDAL_HTTP_COOKIEJAR=/tmp/cookies.txt
# !GDAL_HTTP_COOKIEFILE=/tmp/cookies.txt

# from osgeo import gdal
# gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/cookies.txt')
# gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/cookies.txt')
# gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','EMPTY_DIR') #gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','YES') EMPTY_DIR
# gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF')
# gdal.SetConfigOption('GDAL_HTTP_UNSAFESSL', 'YES')
# gdal.SetConfigOption('GDAL_HTTP_NETRC','True')

# !CPL_VSIL_CURL_USE_HEAD=FALSE # !GDAL_DISABLE_READDIR_ON_OPEN=YES # !GDAL_HTTP_COOKIEJAR=/tmp/cookies.txt # !GDAL_HTTP_COOKIEFILE=/tmp/cookies.txt # from osgeo import gdal # gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/cookies.txt') # gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/cookies.txt') # gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','EMPTY_DIR') #gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','YES') EMPTY_DIR # gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF') # gdal.SetConfigOption('GDAL_HTTP_UNSAFESSL', 'YES') # gdal.SetConfigOption('GDAL_HTTP_NETRC','True')

In [31]:

hls = easysnowdata.remote_sensing.HLS(
    bbox_input=bbox_gdf,
    start_date="2022-07-01",
    end_date="2022-07-31",
    resolution=60,
    add_metadata=False,
    add_platform=False,  # turned off due to bug from https://forum.earthdata.nasa.gov/viewtopic.php?p=20578#p20578, waiting for metadata to be added back
)

#hls.data

hls = easysnowdata.remote_sensing.HLS( bbox_input=bbox_gdf, start_date="2022-07-01", end_date="2022-07-31", resolution=60, add_metadata=False, add_platform=False, # turned off due to bug from https://forum.earthdata.nasa.gov/viewtopic.php?p=20578#p20578, waiting for metadata to be added back ) #hls.data

Data searched. Access the returned seach with the .search_landsat or .search_sentinel attribute.
Data retrieved. Access with the .data attribute. Data CRS: WGS 84 / UTM zone 10N.
Nodata values removed from the data. In doing so, all bands converted to float32. To turn this behavior off, set remove_nodata=False.
Data scaled to reflectance. Access with the .data attribute. To turn this behavior off, set scale_data=False.

In [32]:

hls.get_rgb()
#hls.rgb

hls.get_rgb() #hls.rgb

RGB data retrieved.
Access with the following attributes:
.rgb for raw RGB,
.rgba for RGBA,
.rgb_percentile for percentile RGB,
.rgb_clahe for CLAHE RGB.
You can pass in percentile_kwargs and clahe_kwargs to adjust RGB calculations, check documentation for options.

In [33]:

f,axs = plt.subplots(1,3,figsize=(12,7))

hls.rgba.isel(time=0).plot.imshow(ax=axs[0])
axs[0].set_title('RGBA')
hls.rgb_percentile.isel(time=0).plot.imshow(ax=axs[1])
axs[1].set_title('RGB w/ percentile stretch')
hls.rgb_clahe.isel(time=0).plot.imshow(ax=axs[2])
axs[2].set_title('RGB w/ clahe equalization')

for ax in axs:
    ax.axis('off')
    ax.set_aspect('equal')

f.tight_layout()

f,axs = plt.subplots(1,3,figsize=(12,7)) hls.rgba.isel(time=0).plot.imshow(ax=axs[0]) axs[0].set_title('RGBA') hls.rgb_percentile.isel(time=0).plot.imshow(ax=axs[1]) axs[1].set_title('RGB w/ percentile stretch') hls.rgb_clahe.isel(time=0).plot.imshow(ax=axs[2]) axs[2].set_title('RGB w/ clahe equalization') for ax in axs: ax.axis('off') ax.set_aspect('equal') f.tight_layout()

In [41]:

f = hls.rgb_percentile.plot.imshow(col='time',col_wrap=5,)

# for ax, time, platform in zip(f.axes.flat, hls.rgb.time.values, hls.rgb.platform.values): commented out because of HLS metadata issue that makes platform retrieval impossible, that is being fixed: https://forum.earthdata.nasa.gov/viewtopic.php?p=20578#p20578
#     local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles')
#     ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}\n{platform}')
#     ax.axis('off')
#     ax.set_aspect('equal')

for ax, time in zip(f.axes.flat, hls.rgb_percentile.time.values):
    local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles')
    ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}')
    ax.axis('off')
    ax.set_aspect('equal')

f.fig.subplots_adjust(hspace=0.3)
f.fig.suptitle('Harmonized Landsat Sentinel-2 (HLS) time series\nw/ percentile stretch',fontsize=16,y=1.04)

f = hls.rgb_percentile.plot.imshow(col='time',col_wrap=5,) # for ax, time, platform in zip(f.axes.flat, hls.rgb.time.values, hls.rgb.platform.values): commented out because of HLS metadata issue that makes platform retrieval impossible, that is being fixed: https://forum.earthdata.nasa.gov/viewtopic.php?p=20578#p20578 # local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles') # ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}\n{platform}') # ax.axis('off') # ax.set_aspect('equal') for ax, time in zip(f.axes.flat, hls.rgb_percentile.time.values): local_time = pd.to_datetime(time).tz_localize('UTC').tz_convert('America/Los_Angeles') ax.set_title(f'{local_time.strftime("%B %d, %Y")}\n{local_time.strftime("%I:%M%p")}') ax.axis('off') ax.set_aspect('equal') f.fig.subplots_adjust(hspace=0.3) f.fig.suptitle('Harmonized Landsat Sentinel-2 (HLS) time series\nw/ percentile stretch',fontsize=16,y=1.04)

Out[41]:

Text(0.5, 1.04, 'Harmonized Landsat Sentinel-2 (HLS) time series\nw/ percentile stretch')

MODIS/Terra Snow Cover: MOD10A1 (Daily), MOD10A2 (8-Day), MOD10A1F (Daily cloud-gap-filled)¶

In [42]:

mod10a1 = easysnowdata.remote_sensing.MODIS_snow(bbox_gdf, start_date='2022-07-01', end_date='2022-07-31', data_product='MOD10A1')

mod10a1 = easysnowdata.remote_sensing.MODIS_snow(bbox_gdf, start_date='2022-07-01', end_date='2022-07-31', data_product='MOD10A1')

Data retrieved. Access with the .data attribute.

In [43]:

mod10a1.data

mod10a1.data

Out[43]:

<xarray.Dataset> Size: 4MB
Dimensions:                             (y: 66, x: 167, time: 31)
Coordinates:
  * y                                   (y) float64 528B 5.226e+06 ... 5.196e+06
  * x                                   (x) float64 1kB -9.294e+06 ... -9.217...
    spatial_ref                         int32 4B 0
  * time                                (time) datetime64[ns] 248B 2022-07-01...
Data variables:
    hdf                                 (time, y, x) float32 1MB dask.array<chunksize=(1, 66, 167), meta=np.ndarray>
    NDSI                                (time, y, x) int16 683kB dask.array<chunksize=(1, 66, 167), meta=np.ndarray>
    orbit_pnt                           (time, y, x) uint8 342kB dask.array<chunksize=(1, 66, 167), meta=np.ndarray>
    granule_pnt                         (time, y, x) uint8 342kB dask.array<chunksize=(1, 66, 167), meta=np.ndarray>
    NDSI_Snow_Cover                     (time, y, x) uint8 342kB dask.array<chunksize=(1, 66, 167), meta=np.ndarray>
    Snow_Albedo_Daily_Tile              (time, y, x) uint8 342kB dask.array<chunksize=(1, 66, 167), meta=np.ndarray>
    NDSI_Snow_Cover_Basic_QA            (time, y, x) uint8 342kB dask.array<chunksize=(1, 66, 167), meta=np.ndarray>
    NDSI_Snow_Cover_Algorithm_Flags_QA  (time, y, x) uint8 342kB dask.array<chunksize=(1, 66, 167), meta=np.ndarray>

xarray.Dataset

• Dimensions:
  • y: 66
  • x: 167
  • time: 31
• Coordinates: (4)
  • y
    (y)
    float64
    5.226e+06 5.225e+06 ... 5.196e+06

    units :

    Meter

    resolution :

    -463.31271652750013

    crs :

    PROJCS["unnamed",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",SPHEROID["Custom spheroid",6371007.181,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    array([5225935.786072, 5225472.473355, 5225009.160639, 5224545.847922,
           5224082.535206, 5223619.222489, 5223155.909773, 5222692.597056,
           5222229.28434 , 5221765.971623, 5221302.658907, 5220839.34619 ,
           5220376.033474, 5219912.720757, 5219449.408041, 5218986.095324,
           5218522.782607, 5218059.469891, 5217596.157174, 5217132.844458,
           5216669.531741, 5216206.219025, 5215742.906308, 5215279.593592,
           5214816.280875, 5214352.968159, 5213889.655442, 5213426.342726,
           5212963.030009, 5212499.717293, 5212036.404576, 5211573.09186 ,
           5211109.779143, 5210646.466427, 5210183.15371 , 5209719.840993,
           5209256.528277, 5208793.21556 , 5208329.902844, 5207866.590127,
           5207403.277411, 5206939.964694, 5206476.651978, 5206013.339261,
           5205550.026545, 5205086.713828, 5204623.401112, 5204160.088395,
           5203696.775679, 5203233.462962, 5202770.150246, 5202306.837529,
           5201843.524813, 5201380.212096, 5200916.899379, 5200453.586663,
           5199990.273946, 5199526.96123 , 5199063.648513, 5198600.335797,
           5198137.02308 , 5197673.710364, 5197210.397647, 5196747.084931,
           5196283.772214, 5195820.459498])

  • x
    (x)
    float64
    -9.294e+06 ... -9.217e+06

    units :

    Meter

    resolution :

    463.3127165279165

    crs :

    PROJCS["unnamed",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",SPHEROID["Custom spheroid",6371007.181,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    array([-9294284.749908, -9293821.437192, -9293358.124475, -9292894.811759,
           -9292431.499042, -9291968.186326, -9291504.873609, -9291041.560893,
           -9290578.248176, -9290114.93546 , -9289651.622743, -9289188.310026,
           -9288724.99731 , -9288261.684593, -9287798.371877, -9287335.05916 ,
           -9286871.746444, -9286408.433727, -9285945.121011, -9285481.808294,
           -9285018.495578, -9284555.182861, -9284091.870145, -9283628.557428,
           -9283165.244712, -9282701.931995, -9282238.619279, -9281775.306562,
           -9281311.993845, -9280848.681129, -9280385.368412, -9279922.055696,
           -9279458.742979, -9278995.430263, -9278532.117546, -9278068.80483 ,
           -9277605.492113, -9277142.179397, -9276678.86668 , -9276215.553964,
           -9275752.241247, -9275288.928531, -9274825.615814, -9274362.303098,
           -9273898.990381, -9273435.677665, -9272972.364948, -9272509.052231,
           -9272045.739515, -9271582.426798, -9271119.114082, -9270655.801365,
           -9270192.488649, -9269729.175932, -9269265.863216, -9268802.550499,
           -9268339.237783, -9267875.925066, -9267412.61235 , -9266949.299633,
           -9266485.986917, -9266022.6742  , -9265559.361484, -9265096.048767,
           -9264632.73605 , -9264169.423334, -9263706.110617, -9263242.797901,
           -9262779.485184, -9262316.172468, -9261852.859751, -9261389.547035,
           -9260926.234318, -9260462.921602, -9259999.608885, -9259536.296169,
           -9259072.983452, -9258609.670736, -9258146.358019, -9257683.045303,
           -9257219.732586, -9256756.41987 , -9256293.107153, -9255829.794436,
           -9255366.48172 , -9254903.169003, -9254439.856287, -9253976.54357 ,
           -9253513.230854, -9253049.918137, -9252586.605421, -9252123.292704,
           -9251659.979988, -9251196.667271, -9250733.354555, -9250270.041838,
           -9249806.729122, -9249343.416405, -9248880.103689, -9248416.790972,
           -9247953.478255, -9247490.165539, -9247026.852822, -9246563.540106,
           -9246100.227389, -9245636.914673, -9245173.601956, -9244710.28924 ,
           -9244246.976523, -9243783.663807, -9243320.35109 , -9242857.038374,
           -9242393.725657, -9241930.412941, -9241467.100224, -9241003.787508,
           -9240540.474791, -9240077.162075, -9239613.849358, -9239150.536641,
           -9238687.223925, -9238223.911208, -9237760.598492, -9237297.285775,
           -9236833.973059, -9236370.660342, -9235907.347626, -9235444.034909,
           -9234980.722193, -9234517.409476, -9234054.09676 , -9233590.784043,
           -9233127.471327, -9232664.15861 , -9232200.845894, -9231737.533177,
           -9231274.22046 , -9230810.907744, -9230347.595027, -9229884.282311,
           -9229420.969594, -9228957.656878, -9228494.344161, -9228031.031445,
           -9227567.718728, -9227104.406012, -9226641.093295, -9226177.780579,
           -9225714.467862, -9225251.155146, -9224787.842429, -9224324.529713,
           -9223861.216996, -9223397.904279, -9222934.591563, -9222471.278846,
           -9222007.96613 , -9221544.653413, -9221081.340697, -9220618.02798 ,
           -9220154.715264, -9219691.402547, -9219228.089831, -9218764.777114,
           -9218301.464398, -9217838.151681, -9217374.838965])

  • spatial_ref
    ()
    int32
    0

    spatial_ref :

    PROJCRS["unnamed",BASEGEOGCRS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",ELLIPSOID["Custom spheroid",6371007.181,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["unnamed",METHOD["Sinusoidal"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["Meter",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["Meter",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["northing",north,ORDER[2],LENGTHUNIT["Meter",1]]]

    crs_wkt :

    PROJCRS["unnamed",BASEGEOGCRS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",ELLIPSOID["Custom spheroid",6371007.181,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["unnamed",METHOD["Sinusoidal"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["Meter",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["Meter",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["northing",north,ORDER[2],LENGTHUNIT["Meter",1]]]

    semi_major_axis :

    6371007.181

    semi_minor_axis :

    6371007.181

    inverse_flattening :

    0.0

    reference_ellipsoid_name :

    Custom spheroid

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    Unknown datum based upon the custom spheroid

    horizontal_datum_name :

    Not specified (based on custom spheroid)

    projected_crs_name :

    unnamed

    grid_mapping_name :

    sinusoidal

    longitude_of_projection_origin :

    0.0

    false_easting :

    0.0

    false_northing :

    0.0

    GeoTransform :

    -9294516.406266532838344573974609 463.31271652791650694780401 0 5226167.442430201917886734008789 0 -463.312716527500128904648591

    array(0, dtype=int32)

  • time
    (time)
    datetime64[ns]
    2022-07-01 ... 2022-07-31

    array(['2022-07-01T00:00:00.000000000', '2022-07-02T00:00:00.000000000',
           '2022-07-03T00:00:00.000000000', '2022-07-04T00:00:00.000000000',
           '2022-07-05T00:00:00.000000000', '2022-07-06T00:00:00.000000000',
           '2022-07-07T00:00:00.000000000', '2022-07-08T00:00:00.000000000',
           '2022-07-09T00:00:00.000000000', '2022-07-10T00:00:00.000000000',
           '2022-07-11T00:00:00.000000000', '2022-07-12T00:00:00.000000000',
           '2022-07-13T00:00:00.000000000', '2022-07-14T00:00:00.000000000',
           '2022-07-15T00:00:00.000000000', '2022-07-16T00:00:00.000000000',
           '2022-07-17T00:00:00.000000000', '2022-07-18T00:00:00.000000000',
           '2022-07-19T00:00:00.000000000', '2022-07-20T00:00:00.000000000',
           '2022-07-21T00:00:00.000000000', '2022-07-22T00:00:00.000000000',
           '2022-07-23T00:00:00.000000000', '2022-07-24T00:00:00.000000000',
           '2022-07-25T00:00:00.000000000', '2022-07-26T00:00:00.000000000',
           '2022-07-27T00:00:00.000000000', '2022-07-28T00:00:00.000000000',
           '2022-07-29T00:00:00.000000000', '2022-07-30T00:00:00.000000000',
           '2022-07-31T00:00:00.000000000'], dtype='datetime64[ns]')

• Data variables: (8)
  • hdf
    (time, y, x)
    float32
    dask.array<chunksize=(1, 66, 167), meta=np.ndarray>

    Array Chunk Bytes 1.30 MiB 43.05 kiB Shape (31, 66, 167) (1, 66, 167) Dask graph 31 chunks in 3 graph layers Data type float32 numpy.ndarray

  • NDSI
    (time, y, x)
    int16
    dask.array<chunksize=(1, 66, 167), meta=np.ndarray>

    Array Chunk Bytes 667.35 kiB 21.53 kiB Shape (31, 66, 167) (1, 66, 167) Dask graph 31 chunks in 3 graph layers Data type int16 numpy.ndarray

  • orbit_pnt
    (time, y, x)
    uint8
    dask.array<chunksize=(1, 66, 167), meta=np.ndarray>

    Array Chunk Bytes 333.67 kiB 10.76 kiB Shape (31, 66, 167) (1, 66, 167) Dask graph 31 chunks in 3 graph layers Data type uint8 numpy.ndarray

  • granule_pnt
    (time, y, x)
    uint8
    dask.array<chunksize=(1, 66, 167), meta=np.ndarray>

    Array Chunk Bytes 333.67 kiB 10.76 kiB Shape (31, 66, 167) (1, 66, 167) Dask graph 31 chunks in 3 graph layers Data type uint8 numpy.ndarray

  • NDSI_Snow_Cover
    (time, y, x)
    uint8
    dask.array<chunksize=(1, 66, 167), meta=np.ndarray>

    Array Chunk Bytes 333.67 kiB 10.76 kiB Shape (31, 66, 167) (1, 66, 167) Dask graph 31 chunks in 3 graph layers Data type uint8 numpy.ndarray

  • Snow_Albedo_Daily_Tile
    (time, y, x)
    uint8
    dask.array<chunksize=(1, 66, 167), meta=np.ndarray>

    Array Chunk Bytes 333.67 kiB 10.76 kiB Shape (31, 66, 167) (1, 66, 167) Dask graph 31 chunks in 3 graph layers Data type uint8 numpy.ndarray

  • NDSI_Snow_Cover_Basic_QA
    (time, y, x)
    uint8
    dask.array<chunksize=(1, 66, 167), meta=np.ndarray>

    Array Chunk Bytes 333.67 kiB 10.76 kiB Shape (31, 66, 167) (1, 66, 167) Dask graph 31 chunks in 3 graph layers Data type uint8 numpy.ndarray

  • NDSI_Snow_Cover_Algorithm_Flags_QA
    (time, y, x)
    uint8
    dask.array<chunksize=(1, 66, 167), meta=np.ndarray>

    Array Chunk Bytes 333.67 kiB 10.76 kiB Shape (31, 66, 167) (1, 66, 167) Dask graph 31 chunks in 3 graph layers Data type uint8 numpy.ndarray

• Indexes: (3)
  • y
    PandasIndex

    PandasIndex(Index([ 5225935.786071938,  5225472.473355411,  5225009.160638884,
            5224545.847922356,  5224082.535205828,  5223619.222489301,
            5223155.909772773,  5222692.597056246, 5222229.2843397185,
            5221765.971623191,  5221302.658906664,  5220839.346190136,
            5220376.033473608,  5219912.720757081,  5219449.408040553,
            5218986.095324026,  5218522.782607499,  5218059.469890971,
            5217596.157174444,  5217132.844457916,  5216669.531741388,
            5216206.219024861,  5215742.906308333,  5215279.593591806,
            5214816.280875279,  5214352.968158751,  5213889.655442224,
            5213426.342725696,  5212963.030009168,  5212499.717292641,
            5212036.404576113,  5211573.091859586,  5211109.779143059,
            5210646.466426531,  5210183.153710003,  5209719.840993476,
            5209256.528276948,  5208793.215560421, 5208329.9028438935,
            5207866.590127366,  5207403.277410839,  5206939.964694311,
            5206476.651977783,  5206013.339261256,  5205550.026544728,
            5205086.713828201,  5204623.401111674,  5204160.088395146,
            5203696.775678619,  5203233.462962091,  5202770.150245563,
            5202306.837529036,  5201843.524812508, 5201380.2120959805,
            5200916.899379454,  5200453.586662926,  5199990.273946399,
            5199526.961229871,  5199063.648513343,  5198600.335796816,
            5198137.023080288,  5197673.710363761,  5197210.397647234,
            5196747.084930706,  5196283.772214178,  5195820.459497651],
          dtype='float64', name='y'))

  • x
    PandasIndex

    PandasIndex(Index([-9294284.749908268,  -9293821.43719174, -9293358.124475213,
           -9292894.811758684, -9292431.499042157, -9291968.186325628,
           -9291504.873609101, -9291041.560892573, -9290578.248176046,
           -9290114.935459517,
           ...
           -9221544.653413385, -9221081.340696858,  -9220618.02798033,
           -9220154.715263803, -9219691.402547274, -9219228.089830747,
           -9218764.777114218,  -9218301.46439769, -9217838.151681162,
           -9217374.838964634],
          dtype='float64', name='x', length=167))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2022-07-01', '2022-07-02', '2022-07-03', '2022-07-04',
                   '2022-07-05', '2022-07-06', '2022-07-07', '2022-07-08',
                   '2022-07-09', '2022-07-10', '2022-07-11', '2022-07-12',
                   '2022-07-13', '2022-07-14', '2022-07-15', '2022-07-16',
                   '2022-07-17', '2022-07-18', '2022-07-19', '2022-07-20',
                   '2022-07-21', '2022-07-22', '2022-07-23', '2022-07-24',
                   '2022-07-25', '2022-07-26', '2022-07-27', '2022-07-28',
                   '2022-07-29', '2022-07-30', '2022-07-31'],
                  dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (0)

In [44]:

f= mod10a1.data['NDSI_Snow_Cover'].rio.reproject_match(hls.data, resampling=rio.enums.Resampling.nearest).where(lambda x: x <= 100).plot.imshow(
    col="time", col_wrap=6, vmin=0, vmax=100, cmap='Blues'
)

for ax, time, in zip(f.axs.flat, mod10a1.data.time.values):
    ax.set_title(f'{pd.to_datetime(time).strftime("%B %d, %Y")}')
    ax.axis('off')
    ax.set_aspect('equal')

f.fig.suptitle('MODIS MOD10A1 daily NDSI snow cover',fontsize=16,y=1.02)

f= mod10a1.data['NDSI_Snow_Cover'].rio.reproject_match(hls.data, resampling=rio.enums.Resampling.nearest).where(lambda x: x <= 100).plot.imshow( col="time", col_wrap=6, vmin=0, vmax=100, cmap='Blues' ) for ax, time, in zip(f.axs.flat, mod10a1.data.time.values): ax.set_title(f'{pd.to_datetime(time).strftime("%B %d, %Y")}') ax.axis('off') ax.set_aspect('equal') f.fig.suptitle('MODIS MOD10A1 daily NDSI snow cover',fontsize=16,y=1.02)

Out[44]:

Text(0.5, 1.02, 'MODIS MOD10A1 daily NDSI snow cover')

In [55]:

mod10a2 = easysnowdata.remote_sensing.MODIS_snow(bbox_gdf, start_date='2022-07-01', end_date='2022-07-31', data_product='MOD10A2')

mod10a2 = easysnowdata.remote_sensing.MODIS_snow(bbox_gdf, start_date='2022-07-01', end_date='2022-07-31', data_product='MOD10A2')

Data retrieved. Access with the .data attribute.

In [56]:

mod10a2.get_binary_snow()

mod10a2.get_binary_snow()

Binary snow map calculated. Access with the .binary_snow attribute.

In [57]:

mod10a2.binary_snow.rio.crs

mod10a2.binary_snow.rio.crs

Out[57]:

CRS.from_wkt('PROJCS["unnamed",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",SPHEROID["Custom spheroid",6371007.181,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')

In [58]:

f= mod10a2.binary_snow.rio.reproject_match(hls.data, resampling=rio.enums.Resampling.nearest).plot.imshow(
    col="time", col_wrap=6, cmap='Blues'
)

for ax, time, in zip(f.axs.flat, mod10a2.binary_snow.time.values):
    ax.set_title(f'{pd.to_datetime(time).strftime("%B %d, %Y")}')
    ax.axis('off')
    ax.set_aspect('equal')

f.fig.suptitle('MODIS MOD10A2 8-day snow cover',fontsize=16,y=1.04)

f= mod10a2.binary_snow.rio.reproject_match(hls.data, resampling=rio.enums.Resampling.nearest).plot.imshow( col="time", col_wrap=6, cmap='Blues' ) for ax, time, in zip(f.axs.flat, mod10a2.binary_snow.time.values): ax.set_title(f'{pd.to_datetime(time).strftime("%B %d, %Y")}') ax.axis('off') ax.set_aspect('equal') f.fig.suptitle('MODIS MOD10A2 8-day snow cover',fontsize=16,y=1.04)

Out[58]:

Text(0.5, 1.04, 'MODIS MOD10A2 8-day snow cover')

In [59]:

mod10a1f = easysnowdata.remote_sensing.MODIS_snow(bbox_gdf, start_date='2022-07-01', end_date='2022-07-31', data_product='MOD10A1F')
#NDSI>10 is snow

mod10a1f = easysnowdata.remote_sensing.MODIS_snow(bbox_gdf, start_date='2022-07-01', end_date='2022-07-31', data_product='MOD10A1F') #NDSI>10 is snow

Granules found: 31
 Getting 31 granules, approx download size: 0.03 GB

QUEUEING TASKS | :   0%|          | 0/31 [00:00<?, ?it/s]

File MOD10A1F.A2022182.h09v04.061.2022184060608.hdf already downloaded
File MOD10A1F.A2022183.h09v04.061.2022185043746.hdf already downloaded
File MOD10A1F.A2022184.h09v04.061.2022186055009.hdf already downloaded
File MOD10A1F.A2022185.h09v04.061.2022187234249.hdf already downloaded
File MOD10A1F.A2022186.h09v04.061.2022188071311.hdf already downloaded
File MOD10A1F.A2022187.h09v04.061.2022189044231.hdf already downloaded
File MOD10A1F.A2022188.h09v04.061.2022190054810.hdf already downloaded
File MOD10A1F.A2022189.h09v04.061.2022191062150.hdf already downloaded
File MOD10A1F.A2022190.h09v04.061.2022192033749.hdf already downloaded
File MOD10A1F.A2022191.h09v04.061.2022193061917.hdf already downloaded
File MOD10A1F.A2022192.h09v04.061.2022194154602.hdf already downloaded
File MOD10A1F.A2022193.h09v04.061.2022195213955.hdf already downloaded
File MOD10A1F.A2022194.h09v04.061.2022196065818.hdf already downloaded
File MOD10A1F.A2022195.h09v04.061.2022197050019.hdf already downloaded
File MOD10A1F.A2022196.h09v04.061.2022198043034.hdf already downloaded
File MOD10A1F.A2022197.h09v04.061.2022199044930.hdf already downloaded
File MOD10A1F.A2022198.h09v04.061.2022201021512.hdf already downloaded
File MOD10A1F.A2022199.h09v04.061.2022201063018.hdf already downloaded
File MOD10A1F.A2022200.h09v04.061.2022202070115.hdf already downloaded
File MOD10A1F.A2022202.h09v04.061.2022204074301.hdf already downloaded
File MOD10A1F.A2022201.h09v04.061.2022203064434.hdf already downloaded
File MOD10A1F.A2022204.h09v04.061.2022206074037.hdf already downloaded
File MOD10A1F.A2022205.h09v04.061.2022207080913.hdf already downloaded
File MOD10A1F.A2022206.h09v04.061.2022208051744.hdf already downloaded
File MOD10A1F.A2022207.h09v04.061.2022209055609.hdf already downloaded
File MOD10A1F.A2022203.h09v04.061.2022205045351.hdf already downloaded
File MOD10A1F.A2022208.h09v04.061.2022210065049.hdf already downloaded
File MOD10A1F.A2022209.h09v04.061.2022215090148.hdf already downloaded
File MOD10A1F.A2022210.h09v04.061.2022215090828.hdf already downloaded
File MOD10A1F.A2022211.h09v04.061.2022215091517.hdf already downloaded
File MOD10A1F.A2022212.h09v04.061.2022215092248.hdf already downloaded

PROCESSING TASKS | :   0%|          | 0/31 [00:00<?, ?it/s]

COLLECTING RESULTS | :   0%|          | 0/31 [00:00<?, ?it/s]

Data retrieved. Access with the .data attribute.

In [60]:

f= mod10a1f.data.rio.reproject_match(hls.data, resampling=rio.enums.Resampling.nearest).where(lambda x: x <= 100).plot.imshow(
    col="time", col_wrap=6, vmin=0, vmax=100, cmap='Blues'
)

for ax, time, in zip(f.axs.flat, mod10a1f.data.time.values):
    ax.set_title(f'{pd.to_datetime(time).strftime("%B %d, %Y")}')
    ax.axis('off')
    ax.set_aspect('equal')

f.fig.suptitle('MODIS MOD10A1F cloud-gap-filled NDSI snow cover',fontsize=16,y=1.02)

f= mod10a1f.data.rio.reproject_match(hls.data, resampling=rio.enums.Resampling.nearest).where(lambda x: x <= 100).plot.imshow( col="time", col_wrap=6, vmin=0, vmax=100, cmap='Blues' ) for ax, time, in zip(f.axs.flat, mod10a1f.data.time.values): ax.set_title(f'{pd.to_datetime(time).strftime("%B %d, %Y")}') ax.axis('off') ax.set_aspect('equal') f.fig.suptitle('MODIS MOD10A1F cloud-gap-filled NDSI snow cover',fontsize=16,y=1.02)

Out[60]:

Text(0.5, 1.02, 'MODIS MOD10A1F cloud-gap-filled NDSI snow cover')

PALSAR-2¶

In [43]:

# list(bbox_gdf.total_bounds)

# !pip install --upgrade xee

# import xee

# !earthengine authenticate --quiet

# import ee

# ee.Initialize(opt_url="https://earthengine-highvolume.googleapis.com")

# bbox_ee = ee.Geometry.Rectangle(153.0, -43.0, 154.0, -42.0)
# bbox_ee = ee.Geometry.Rectangle(*list(bbox_gdf.total_bounds))

# bbox_ee = ee.Geometry.BBox(*list(bbox_gdf.total_bounds))

# bbox_ee.Polygon()

# ic.first().select(0).projection()

# ic = ee.ImageCollection("JAXA/ALOS/PALSAR-2/Level2_2/ScanSAR").filterDate(
#     "2022-07-01", "2022-07-31"
# )
# ds = xr.open_dataset(
#     ic,
#     geometry=bbox_ee,
#     projection=ic.first().select(0).projection(),
#     crs="EPSG:4326",
#     engine="ee",
#     scale=0.01,
# )

# ds

# ds["HH"].isel(time=1).plot()

# def dn_to_db(dn):
#     db = 10 * np.log10(dn**2) - 83
#     return db

# ds["HH"].map_blocks(dn_to_db)
# ds["HH"].isel(time=0)
# ds.time.values[0:10]

# list(bbox_gdf.total_bounds) # !pip install --upgrade xee # import xee # !earthengine authenticate --quiet # import ee # ee.Initialize(opt_url="https://earthengine-highvolume.googleapis.com") # bbox_ee = ee.Geometry.Rectangle(153.0, -43.0, 154.0, -42.0) # bbox_ee = ee.Geometry.Rectangle(*list(bbox_gdf.total_bounds)) # bbox_ee = ee.Geometry.BBox(*list(bbox_gdf.total_bounds)) # bbox_ee.Polygon() # ic.first().select(0).projection() # ic = ee.ImageCollection("JAXA/ALOS/PALSAR-2/Level2_2/ScanSAR").filterDate( # "2022-07-01", "2022-07-31" # ) # ds = xr.open_dataset( # ic, # geometry=bbox_ee, # projection=ic.first().select(0).projection(), # crs="EPSG:4326", # engine="ee", # scale=0.01, # ) # ds # ds["HH"].isel(time=1).plot() # def dn_to_db(dn): # db = 10 * np.log10(dn**2) - 83 # return db # ds["HH"].map_blocks(dn_to_db) # ds["HH"].isel(time=0) # ds.time.values[0:10]