diff --git a/azure-landsat/landsat-hls-azure.py b/azure-landsat/landsat-hls-azure.py index 3322ce4..088022c 100755 --- a/azure-landsat/landsat-hls-azure.py +++ b/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,11 +76,10 @@ 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 @@ -87,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' @@ -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)