diff --git a/azure-landsat/landsat-c2.py b/azure-landsat/landsat-c2.py index e62fac4..98b4489 100755 --- a/azure-landsat/landsat-c2.py +++ b/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) diff --git a/azure-landsat/setup.sh b/azure-landsat/setup.sh index 2af4a10..483c963 100755 --- a/azure-landsat/setup.sh +++ b/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