Option for Snappy ProductIO module to read the bands of tiff image-CodePudding

In the

But I failed to install this module to my Python 3.8 environment. Instead, I tried using a different package to read the bands of the 31-band tiff image.

from osgeo import gdal
from PIL import Image
import numpy as np
import matplotlib as mtp
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import earthpy.plot as ep
import rasterio
from rasterio.plot import reshape_as_raster, reshape_as_image

%matplotlib inline

pd.options.display.max_colwidth = 89

#setting the path for image
S1_S2_stack = 'S1_S2_stack.tif'

#path to training and validation data
training_points = 'testing.shp'
validation_points = 'training.shp'

colors = dict ((
    (0, (0,76,153,255)),  #wheat
    (1, (0,153,0,255)),   #corn
    (2, (255,0,0,255)),   #other
    (3, (255,153,51,255)),
    (4, (255,255,0,255))
    
))

for k in colors:
    v = colors [k]
    _v = [_v / 255.0 for _v in v]
    colors[k] = _v
    
index_colors = [colors[key] if key in colors else (1,1,1,0) for key in range (0,5)]
cmap = plt.matplotlib.colors.ListedColormap(index_colors, 'Classification', 5)

src = rasterio.open(S1_S2_stack)
src1 = src.read(S1_S2_stack)
bands = list (src1.tags())

And when I run the last section it throws me an error which says:

IndexError: band index S out of range (not in (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31))

So I'll appreciate any other suggestions.

CodePudding user response：

The file name is not a band index.

S1_S2_stack = 'S1_S2_stack.tif'
src = rasterio.open(S1_S2_stack)
src1 = src.read(S1_S2_stack)      # Error: sending a name instead of an index

Here's what we know about the read method:

def read(indexes=None, ...)
"""
indexes : int or list, optional
    If `indexes` is a list, the result is a 3D array, but is
    a 2D array if it is a band index number.
"""

It reads the bands by their indexes, which are the first parameter and should be either an integer or a list of integers. If None it returns all the bands. Therefore the code should be something like

src1 = src.read()            # Read all bands together

index = 0                    # set the band to read
assert 1 <= index <= src.count    # make sure the index does not exceed the number of bands
src1 = src.read(index)       # read the band

indexes = [1,2,3]
src1 = src.read(indexes)     # read some of the bands

CodePudding user response：

I tried your suggestion @Vitalizzare. I think I need to share my code completely to find the answer.

from osgeo import gdal
from PIL import Image
import numpy as np
import matplotlib as mtp
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import earthpy.plot as ep
import rasterio
from rasterio.plot import reshape_as_raster, reshape_as_image

get_ipython().run_line_magic('matplotlib', 'inline')

pd.options.display.max_colwidth = 89


# In[2]:


#setting the path for image
S1_S2_stack = 'S1_S2_stack.tif'

#path to training and validation data
training_points = 'testing.shp'
validation_points = 'training.shp'


# In[3]:


colors = dict ((
    (0, (0,76,153,255)),  #wheat
    (1, (0,153,0,255)),   #corn
    (2, (255,0,0,255)),   #other
    (3, (255,153,51,255)),
    (4, (255,255,0,255))
    
))


# In[4]:


for k in colors:
    v = colors [k]
    _v = [_v / 255.0 for _v in v]
    colors[k] = _v
    
index_colors = [colors[key] if key in colors else (1,1,1,0) for key in range (0,5)]
cmap = plt.matplotlib.colors.ListedColormap(index_colors, 'Classification', 5)


# In[5]:


src = rasterio.open(S1_S2_stack)
bands = src.read()


# In[6]:


stack =src. read()
print(stack.shape)

fig, (ax1, ax2) = plt.subplots(1,2,figsize= (20,10))
ax1 = ep.plot_rgb(arr = stack, rgb =(3, 2, 1), stretch=True, ax = ax1, title = "RGB Composite - Sentinel-2")
stack_s1 =np.stack ((stack[28],stack[29],stack[29]/stack[28]))
ax2 = ep.plot_rgb(arr = stack_s1, rgb = (1,0,2), stretch=True, ax = ax2, title= "RGB Composite - Sentinel-1 (VV, VH, VV/VH)")
plt.tight_layout()


# In[7]:


img = src.read()
profile = src.profile
with rasterio.io.MemoryFile () as memfile:
    with memfile.open(**profile) as dst:
        for i in range(0, src.count):
            dst.write(img[i], i 1)
    dataset = memfile.open()

Works great from here. But this one gave an error:

# In[8]:


#read points from shapefile
train_pts = gpd.read_file (training_points)
train_pts = train_pts[[ 'CID','class', 'POINT_X','POINT_Y']] #attribute fields in shapefile
train_pts.index = range(len(train_pts))
coords = [(x,y) for x, y in zip(train_pts.POINT_X, train_pts.POINT_Y)] #create list of point coordinates

#sample each band of raster dataset at each point in the coordinate list
train_pts ['Raster Value'] = [x for x in dataset.sample(coords)] #all band values saved as a list in the Raster Value column
#Unpack the raster value column to separate column for each band - band names retrieved from snappy in the video but I was looking for an option
train_pts[bands] = pd.DataFrame(train_pts['Raster Value'].tolist(), index = train_pts.index)
train_pts = train_pts.drop(['Raster Value'], axis=1) #drop raster value column
#change the values for last three classes 
train_pts['CID'] = train_pts['CID'].replace([7,8,15],[5,6,7])
train_pts.to.csv('train_pts.csv') #save as csv
train_pts.head () #see first column

"ValueError: Array conditional must be same shape as self"