update setup.sh

azure-landsat/landsat-c2.py

@@ -1,7 +1,6 @@
 #!/usr/bin/env python3
 # reference - https://planetarycomputer.microsoft.com/docs/concepts/sas/
 # source - https://github.com/microsoft/AIforEarthDataSets/blob/main/data/landsat-7.ipynb
-# get key curl https://planetarycomputer.microsoft.com/api/sas/v1/token/landsateuwest/landsat-c2 | cut -d '"' -f 8 > file
 
 import os
 import datetime
@@ -155,40 +154,6 @@
 for band_name in rgb_bands:
     rgb_urls.append([s for s in azure_urls if band_name + '.TIF' in s][0])
 
-thumbnail_data = []
-
-# url = rgb_urls[0]
-for url in rgb_urls:
-
-#    print(url)
-
-    # From:
-    #
-    # https://automating-gis-processes.github.io/CSC/notebooks/L5/read-cogs.html
-    with rasterio.open(url) as raster:
-
-        # List of overviews from biggest to smallest
-        oviews = raster.overviews(1)
-
-        # Retrieve a small-ish thumbnail (six are available for Landsat files)
-        decimation_level = oviews[thumbnail_index]
-        h = int(raster.height/decimation_level)
-        w = int(raster.width/decimation_level)
-
-        thumbnail_channel = raster.read(1, out_shape=(1, h, w)) / composite_norm_value
-        thumbnail_data.append(thumbnail_channel)
-
-rgb1 = np.dstack((thumbnail_data[0],thumbnail_data[1],thumbnail_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.imshow(rgb1);
-#plt.show()
-#plt.imsave('test.png', rgb)
-
 image_data = []
 
 for fn in rgb_urls:
@@ -200,50 +165,12 @@
         band_array = band_array / composite_norm_value
         image_data.append(band_array)
 
-rgb2 = np.dstack((image_data[0],image_data[1],image_data[2]))
-np.clip(rgb2,0,1,rgb2)
+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(rgb2);
+#plt.imshow(rgb);
 #plt.show()
-
-image_data = []
-
-for fn in rgb_urls:
-
-    with rasterio.open(fn,'r') as src:
-
-        # Choose the target image size
-        xsize = (target_size[0] * oversample)
-        ysize = (target_size[1] * oversample)
-
-        # Find the transformation from pixels to lat/lon
-        transformer = Transformer.from_crs('EPSG:4326', src.crs, always_xy=True)
-        xx, yy = transformer.transform(query_lon, query_lat)
-        py, px = src.index(xx,yy)
-        xoff = px - xsize // 2 
-        yoff = py - ysize // 2 
-        window = Window(xoff,yoff,xsize,ysize)
-
-        band_array = src.read(1, window=window)
-        src.close()
-        band_array = band_array / composite_norm_value
-        image_data.append(band_array)
-
-    # ...with rasterio.open()
-
-# ...for each file
-
-rgb3 = np.dstack((image_data[0],image_data[1],image_data[2]))
-np.clip(rgb3,0,1,rgb3)
-
-w = target_size[0]
-h = target_size[1]
-
-dpi = 50; 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(rgb3);
-plt.show();
+plt.imsave('Shizuoka-Japan.png', rgb)

azure-landsat/setup.sh

@@ -16,6 +16,5 @@ if [ ! -d ./tokens ]; then
   mkdir -p ./tokens;
 fi
 
-key=$(curl -s https://planetarycomputer.microsoft.com/api/sas/v1/token/landsateuwest/landsat-c2 | cut -d '"' -f 8 >  tokens/file)
-echo $key
-
+key=$(curl -s https://planetarycomputer.microsoft.com/api/sas/v1/token/landsateuwest/landsat-c2 | cut -d '"' -f 8 > tokens/landsat_ro_sas.txt)
+#echo $key