From a1d8b54bd06faa2f54be6613ffce3aac4a08bd32 Mon Sep 17 00:00:00 2001
From: John Hicks <jhicks@internet2.edu>
Date: Mon, 28 Feb 2022 15:50:17 -0500
Subject: [PATCH 1/2] cleaned up landsat analysis using VScode
---
azure-landsat/landsat-hls-azure.py | 15 +++++++--------
1 file changed, 7 insertions(+), 8 deletions(-)
diff --git a/azure-landsat/landsat-hls-azure.py b/azure-landsat/landsat-hls-azure.py
index 3322ce4..70022ac 100755
--- a/azure-landsat/landsat-hls-azure.py
+++ b/azure-landsat/landsat-hls-azure.py
@@ -1,6 +1,4 @@
#!/usr/bin/env python3
-# need to install - pip install azure-storage-blob
-# https://nbviewer.org/github/microsoft/AIforEarthDataSets/blob/main/data/hls.ipynb
# Environment setup
import requests
@@ -75,8 +73,8 @@ def lat_lon_to_hls_tile_id(lat,lon):
# Build a tile path from components, including lat/lon lookup
# Near Bear Lake, UT
-# lat = 41.89655047211277; lon = -111.4132464403312
-lat = 39.768402; lon = -86.158066
+lat = 41.89655047211277; lon = -111.4132464403312
+# lat = 39.768402; lon = -86.158066
tile_id = lat_lon_to_hls_tile_id(lat,lon)
print('Using Tile ID {}'.format(tile_id))
@@ -120,8 +118,8 @@ def lat_lon_to_hls_tile_id(lat,lon):
dpi = 100; fig = plt.figure(frameon=False,figsize=(w/dpi,h/dpi),dpi=dpi)
ax = plt.Axes(fig,[0., 0., 1., 1.]); ax.set_axis_off(); fig.add_axes(ax)
-plt.imshow(rgb);
-plt.imsave('test.png', rgb)
+# plt.imshow(rgb);
+# plt.imsave('test.png', rgb)
i_daynum = int(daynum)
product = 'L30'
@@ -165,5 +163,6 @@ def lat_lon_to_hls_tile_id(lat,lon):
dpi = 100; fig = plt.figure(frameon=False,figsize=(w/dpi,h/dpi),dpi=dpi)
ax = plt.Axes(fig,[0., 0., 1., 1.]); ax.set_axis_off(); fig.add_axes(ax)
-plt.imshow(rgb);
-plt.show()
+# plt.imshow(rgb);
+# plt.show()
+plt.imsave('test.png', rgb)
From c5659b1a9715034962b06ae7af842c781d5157ce Mon Sep 17 00:00:00 2001
From: John Hicks <jhicks@internet2.edu>
Date: Mon, 7 Mar 2022 17:09:06 -0500
Subject: [PATCH 2/2] updated landsat code
---
azure-landsat/landsat-hls-azure.py | 64 +++++++++++++++---------------
1 file changed, 33 insertions(+), 31 deletions(-)
diff --git a/azure-landsat/landsat-hls-azure.py b/azure-landsat/landsat-hls-azure.py
index 70022ac..088022c 100755
--- a/azure-landsat/landsat-hls-azure.py
+++ b/azure-landsat/landsat-hls-azure.py
@@ -1,4 +1,7 @@
#!/usr/bin/env python3
+# The python scripe originally came from the following jupyter nodebook
+# https://nbviewer.org/github/microsoft/AIforEarthDataSets/blob/main/data/hls.ipynb
+#
# Environment setup
import requests
@@ -31,7 +34,7 @@
# Tile ID, Xstart, Ystart, UZ, EPSG, MinLon, MaxLon, MinLon, MaxLon
s = requests.get(hls_tile_extents_url).content
hls_tile_extents = pd.read_csv(io.StringIO(s.decode('utf-8')),delimiter=r'\s+')
-print('Read tile extents for {} tiles'.format(len(hls_tile_extents)))
+#print('Read tile extents for {} tiles'.format(len(hls_tile_extents)))
hls_container_client = ContainerClient(account_url=storage_account_url,
container_name=container_name,
@@ -74,10 +77,9 @@ def lat_lon_to_hls_tile_id(lat,lon):
# Build a tile path from components, including lat/lon lookup
# Near Bear Lake, UT
lat = 41.89655047211277; lon = -111.4132464403312
-# lat = 39.768402; lon = -86.158066
tile_id = lat_lon_to_hls_tile_id(lat,lon)
-print('Using Tile ID {}'.format(tile_id))
+#print('Using Tile ID {}'.format(tile_id))
year = '2019'
daynum = '001' # 1-indexed day-of-year
@@ -85,41 +87,41 @@ def lat_lon_to_hls_tile_id(lat,lon):
sband = '_01'
version = 'v1.4' # Currently always v1.4
-blob_name = product + '/HLS.' + product + '.T' + tile_id + '.' + year + daynum + '.' + version \
- + sband + '.tif'
-sentinel_url = hls_blob_root + '/' + blob_name
-print('Constructed tile URL:\n{}'.format(sentinel_url))
+#blob_name = product + '/HLS.' + product + '.T' + tile_id + '.' + year + daynum + '.' + version \
+# + sband + '.tif'
+#sentinel_url = hls_blob_root + '/' + blob_name
+#print('Constructed tile URL:\n{}'.format(sentinel_url))
# Display an RGB composite of the Sentinel-2 image we just chose
# Bands 2, 3, and 4 are B, G, and R in Sentinel-2 HLS images
composite_bands = [4,3,2]
-urls = [sentinel_url.replace(sband,'_' + str(b).zfill(2)) for b in composite_bands]
-print('Rendering URLs:')
-for s in urls:
- print(s)
+#urls = [sentinel_url.replace(sband,'_' + str(b).zfill(2)) for b in composite_bands]
+#print('Rendering URLs:')
+#for s in urls:
+# print(s)
-image_data = []
+#image_data = []
composite_downsample_factor = 5
composite_norm_value = 15000
-for fn in urls:
- with rasterio.open(fn,'r') as raster:
- h = int(raster.height/composite_downsample_factor)
- w = int(raster.width/composite_downsample_factor)
- band_array = raster.read(1, out_shape=(1, h, w))
- raster.close()
- band_array = band_array / composite_norm_value
- image_data.append(band_array)
+#for fn in urls:
+# with rasterio.open(fn,'r') as raster:
+# h = int(raster.height/composite_downsample_factor)
+# w = int(raster.width/composite_downsample_factor)
+# band_array = raster.read(1, out_shape=(1, h, w))
+# raster.close()
+# band_array = band_array / composite_norm_value
+# image_data.append(band_array)
-rgb = np.dstack((image_data[0],image_data[1],image_data[2]))
-np.clip(rgb,0,1,rgb)
+#rgb = np.dstack((image_data[0],image_data[1],image_data[2]))
+#np.clip(rgb,0,1,rgb)
-dpi = 100; fig = plt.figure(frameon=False,figsize=(w/dpi,h/dpi),dpi=dpi)
-ax = plt.Axes(fig,[0., 0., 1., 1.]); ax.set_axis_off(); fig.add_axes(ax)
-# plt.imshow(rgb);
-# plt.imsave('test.png', rgb)
+#dpi = 100; fig = plt.figure(frameon=False,figsize=(w/dpi,h/dpi),dpi=dpi)
+#ax = plt.Axes(fig,[0., 0., 1., 1.]); ax.set_axis_off(); fig.add_axes(ax)
+#plt.imshow(rgb);
+#plt.imsave('test01.png', rgb)
i_daynum = int(daynum)
product = 'L30'
@@ -158,11 +160,11 @@ def lat_lon_to_hls_tile_id(lat,lon):
band_array = band_array / composite_norm_value
image_data.append(band_array)
-rgb = np.dstack((image_data[0],image_data[1],image_data[2]))
-np.clip(rgb,0,1,rgb)
+rgb1 = np.dstack((image_data[0],image_data[1],image_data[2]))
+np.clip(rgb1,0,1,rgb1)
dpi = 100; fig = plt.figure(frameon=False,figsize=(w/dpi,h/dpi),dpi=dpi)
ax = plt.Axes(fig,[0., 0., 1., 1.]); ax.set_axis_off(); fig.add_axes(ax)
-# plt.imshow(rgb);
-# plt.show()
-plt.imsave('test.png', rgb)
+#plt.imshow(rgb1);
+#plt.show()
+plt.imsave('test.png', rgb1)