Azure

azure-landsat/landsat-hls-azure.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
-# need to install - pip install azure-storage-blob
+# 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
@@ -33,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,
@@ -75,53 +76,52 @@ 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 
 
 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
 product = 'S30'   # 'S30' for Sentinel, 'L30' for Landsat
 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'
@@ -160,10 +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.imshow(rgb1);
+#plt.show()
+plt.imsave('test.png', rgb1)