<jupyter_start><jupyter_text>Code adapted from https://github.com/patrickcgray/open-geo-tutorial<jupyter_code>from IPython.display import Audio, display
from timeit import default_timer as timer
start = timer()
def alert():
  global start
  end= timer()
  print("\n\nDuration in minutes : " +  str((end - start)/60) + " minutes." )
  display(Audio(url='https://sound.peal.io/ps/audios/000/000/537/original/woo_vu_luvub_dub_dub.wav', autoplay=True))<jupyter_output><empty_output><jupyter_text>Preparing Dataset  Cutting images for faster process<jupyter_code># 0 to  10980 for this img_data
#training
xmin = 6500
xmax = 7000
ymin = 6600
ymax = 7100
#test
xmin_test = 6000
xmax_test = 6500
ymin_test = 6800
ymax_test = 7300
import rasterio
from rasterio.mask import mask
import geopandas as gp
import numpy as np
import shapely
from shapely import geometry
from shapely.geometry import shape, Point, LineString, Polygon , mapping
import matplotlib.pyplot as plt
from rasterio.plot import show
import pyproj
from pyproj import CRS
import fiona
import os # we need os to do some basic file operations
import re
regexp = re.compile(r'_B0[2348]')
sentinal_20m = "./data/IMG_DATA_1"
# find every file in the sentinal_fp directory
sentinal_band_paths = [os.path.join(sentinal_20m, f) for f in os.listdir(sentinal_20m) if regexp.search(f) ]
sentinal_band_paths.sort()
sentinal_band_paths
for i in sentinal_band_paths : 
        src = rasterio.open(i,driver='JP2OpenJPEG')
        meta = src.meta
        print( i + " ::: ")
        print(meta)

# create a products directory within the data dir which won't be uploaded to Github
img_dir = './data/'

# check to see if the dir it exists, if not, create it
if not os.path.exists(img_dir):
    os.makedirs(img_dir)

# filepath for image we're writing out
img_fp = img_dir + 'sentinel_bands_img_data_1.tif'<jupyter_output><empty_output><jupyter_text>Only run this the first time - its slow<jupyter_code>"""
# Just run this the first time 
meta = meta
# Update metadata to reflect the number of layers
meta.update(
    count = len(sentinal_band_paths)
)
# Read each layer and write it to stack
with rasterio.open(img_fp, 'w', **meta) as dst:
    for id, layer in enumerate(sentinal_band_paths, start=1):
        with rasterio.open(layer) as src1:
            dst.write_band(id, src1.read(1))
"""
full_dataset = rasterio.open(img_fp)
img_rows, img_cols = full_dataset.shape
img_bands = full_dataset.count
print(full_dataset.shape) # dimensions
print(full_dataset.count) # bands
dataset=full_dataset
# what is the name of this image
img_name = dataset.name
print('Image filename: {n}\n'.format(n=img_name))

# How many bands does this image have?
num_bands = dataset.count
print('Number of bands in image: {n}\n'.format(n=num_bands))

# How many rows and columns?
rows, cols = dataset.shape
print('Image size is: {r} rows x {c} columns\n'.format(r=rows, c=cols))

# Does the raster have a description or metadata?
desc = dataset.descriptions
metadata = dataset.meta

print('Raster description: {desc}\n'.format(desc=desc))

# What driver was used to open the raster?
driver = dataset.driver
print('Raster driver: {d}\n'.format(d=driver))

# What is the raster's projection?
proj = dataset.crs
if (not proj) :
    proj = CRS.from_epsg(4326)
print('Image projection:')
print(proj, '\n')

# What is the raster's "geo-transform"
gt = dataset.transform

print('Image geo-transform:\n{gt}\n'.format(gt=gt))

print('All raster metadata:')
print(metadata)
print('\n')
from rasterio.plot import adjust_band
from rasterio.plot import reshape_as_raster, reshape_as_image


cutted = dataset.read([3,2,1])[:,xmin: xmax,ymin : ymax]
rgb_norm = adjust_band(cutted) # normalize bands to range between 1.0 to 0.0
show(rgb_norm)

cutted = dataset.read([3,2,1])[:,xmin_test: xmax_test,ymin_test : ymax_test]
rgb_norm = adjust_band(cutted) # normalize bands to range between 1.0 to 0.0

show(rgb_norm)
import matplotlib.pyplot as plt
from rasterio.plot import show
"""
clipped_img = adjust_band(cutted)
print(clipped_img.shape)
fig, ax = plt.subplots(figsize=(10,7))
show(clipped_img[:,  : , : ], ax=ax, transform=full_dataset.transform) # add the transform arg to get it in lat long coords
"""<jupyter_output><empty_output><jupyter_text>Indices<jupyter_code>with rasterio.open(img_fp) as src:
    green_band = src.read(2)[ xmin : xmax ,  ymin : ymax ] # band 3
    red_band = src.read(3)[ xmin : xmax ,  ymin : ymax ]   # band 4
    nir_band = src.read(4)[ xmin : xmax ,  ymin : ymax ]   # band 8
    img = src.read()[:,  xmin : xmax ,  ymin : ymax ]
    
ndwi = (green_band.astype(float) - nir_band.astype(float)) / (green_band.astype(float) + nir_band.astype(float))
ndvi = (nir_band.astype(float) - red_band.astype(float)) / (red_band.astype(float) + nir_band.astype(float))
def color_stretch(image, index):
    colors = image[:, :, index].astype(np.float64)
    for b in range(colors.shape[2]):
        colors[:, :, b] = rasterio.plot.adjust_band(colors[:, :, b])
    return colors
"""
fig, axs = plt.subplots(2,2,figsize=(15,7))

reshaped_img =reshape_as_image( img )
img_stretched = color_stretch(reshaped_img, [2, 1, 0])
axs[0,0].set_title("Stretch")
axs[0,0].imshow(img_stretched)

axs[0,1].set_title("True Color")
axs[0,1].imshow(reshape_as_image(rgb_norm))

nwdi_plot = axs[1,0].imshow(ndwi, cmap="RdYlBu")
axs[1,0].set_title("NDWI")
fig.colorbar(nwdi_plot, ax=axs[1,0])

ndvi_plot = axs[1,1].imshow(ndvi, cmap="RdYlGn")
axs[1,1].set_title("NDVI")
fig.colorbar(ndvi_plot, ax=axs[1,1])

plt.show()
"""
src.close()<jupyter_output><empty_output><jupyter_text>Intercepting label data with cutting image<jupyter_code># Open the dataset from the file
shapefile = gp.read_file('./data/shapefiles/shapefile_compressed/shp_2012_IMG_DATA_1.shp')
shapefile = shapefile.filter(['Legenda','geometry','CLC2012'])
shapefile.crs = {'init': 'epsg:4326'}
shapefile
shapefile = shapefile.to_crs({'init': 'epsg:32629'})
shapefile.bounds<jupyter_output><empty_output><jupyter_text>Img_data_1 cutted<jupyter_code>from pyproj import transform
from pyproj import Proj
def cut_bounds( xmin , xmax ,  ymin , ymax ):
    # this will get our four corner points
    xmin = xmin 
    ymin = ymin
    xmax = xmax
    ymax = ymax
    raster_gps_points = full_dataset.transform * (xmin, ymin),full_dataset.transform * (xmax, ymin),full_dataset.transform * (xmax, ymax), full_dataset.transform * (xmin, ymax),
    # Project all longitudes, latitudes using the pyproj package
    p1 = Proj(dataset.crs) # our current crs
    print(p1)
    print('raster bounds in current crs :\n', raster_gps_points)
    return raster_gps_points
raster_gps_points = cut_bounds( xmin , xmax ,  ymin , ymax )
raster_gps_points_test = cut_bounds( xmin_test , xmax_test ,  ymin_test , ymax_test )
def shapefile_from_rasterbounds(raster_gps_points,shapefile) :
    polygon = Polygon( list(raster_gps_points))
    print("Polygon raster gps points",polygon)
    imageBounds = gp.GeoDataFrame(crs=CRS.from_epsg(32629))
    imageBounds.geometry = [polygon]
    imageBounds = imageBounds.to_crs(fiona.crs.from_epsg(4326))
    print("Imaged cutted bounds",imageBounds.bounds)
    shapefile = shapefile.to_crs(fiona.crs.from_epsg(4326))
    print("Shapefile bounds",shapefile.bounds)
    intersected = gp.sjoin( shapefile, imageBounds, how='inner', op='intersects', lsuffix='left', rsuffix='right')
    intersected = intersected.reset_index()
    print("Intersected df",intersected)
    print("Intersected crs",intersected.crs)
    print("Intersected bounds",intersected.bounds)
    return intersected
intersected = shapefile_from_rasterbounds(raster_gps_points,shapefile)
intersected_test = shapefile_from_rasterbounds(raster_gps_points_test,shapefile)
import re
def class_from_CLCcode(clc) :
    mega_classes = {
        '^1.*' : 'Territórios artificializados',
        '^2.*' : 'Agricultura',
        '^31.*' : 'Floresta',
        '^3.4' : 'Floresta',
        '^32[123]' : 'Vegetação natural',
        '^33[123]' : 'Espaços descobertos ou com vegetação esparsa',
        '^41.*' : 'Zonas húmidas',
        '^42[13]' : 'Zonas húmidas',
        '^5.*' : 'Corpos de água',
        '^422' : 'Corpos de água'
    }
    for mega_class_key in mega_classes :
        searched =  re.search(mega_class_key, str(clc))
        if searched is not None:
            return mega_classes[mega_class_key]
    return 'Not Defined'

x = intersected['CLC2012'][0] 
print (x)
class_from_CLCcode (x) 
for i in intersected.index:
    clc = intersected.at[i, 'CLC2012']
    new_class = class_from_CLCcode(clc)
    intersected.at[i, 'Mega_Legenda'] = new_class
intersected.head()
for i in intersected_test.index:
    clc = intersected_test.at[i, 'CLC2012']
    new_class = class_from_CLCcode(clc)
    intersected_test.at[i, 'Mega_Legenda'] = new_class
intersected_test.head()
import matplotlib.pyplot as plt
unique, counts = np.unique(intersected['Mega_Legenda'], return_counts=True)
plt.bar(unique, counts, 1)
plt.title('Class Frequency')
plt.xlabel('Class')
plt.ylabel('Frequency')
plt.show()

print('Number of different classes :'  + str(len(unique)))
print(unique)<jupyter_output><empty_output><jupyter_text>Setting up label data with bands<jupyter_code>def setup_shafile_to_data(shapefile,full_dataset) :
    shapefile = shapefile.filter(['Legenda','Mega_Legenda', 'geometry','CLC2012'])
    print("Data crs: " + str(full_dataset.crs))
    print("Before crs: " + str(shapefile.crs))
    if shapefile.crs != {'init': 'epsg:32629'} :
        shapefile = shapefile.to_crs({'init': 'epsg:32629'})
    print("After crs: " + str(shapefile.crs))
    return shapefile
shapefile = setup_shafile_to_data(intersected,full_dataset)
shapefile_test = setup_shafile_to_data(intersected_test,full_dataset)
shapefile.head()<jupyter_output><empty_output><jupyter_text>Choose target column = Legenda or Mega_Legenda<jupyter_code>Class = 'Mega_Legenda'
#Class = 'Legenda' 
unique= np.unique(shapefile[Class])
print('List of Land Cover Classes:')
unique
unique= np.unique(shapefile_test[Class])
print('List of Land Cover Classes for test:')
unique
#shapefile.geometry = shapefile.geometry.map(lambda polygon: shapely.ops.transform(lambda x, y: (y, x), polygon))
shapefile.bounds
print('raster bounds in current crs :\n', raster_gps_points)
# this generates a list of shapely geometries
geoms = shapefile.geometry.values 
geoms_test = shapefile_test.geometry.values 

# let's grab a single shapely geometry to check
geometry = geoms[0] 
print(type(geometry))
print(geometry)

# transform to GeoJSON format
from shapely.geometry import mapping
feature = [mapping(geometry)] # can also do this using polygon.__geo_interface__
print(type(feature))
print(feature)
out_image, out_transform = rasterio.mask.mask(full_dataset, feature,crop=True)
out_image.shape
out_transform
def setup_features_and_target(src,geoms,shapefile) : 
    X = np.array([], dtype=np.int8).reshape(0,4) # pixels for training
    y = np.array([], dtype=np.string_) # labels for training
    # extract the raster values within the polygon 
    band_count = src.count
    for index, geom in enumerate(geoms):
        feature = [mapping(geom)]

        # the mask function returns an array of the raster pixels within this feature
        out_image, out_transform = rasterio.mask.mask(src, feature,crop=True)
        # eliminate all the pixels with 0 values for all 8 bands - AKA not actually part of the shapefile
        out_image_trimmed = out_image[:,~np.all(out_image == 0, axis=0)]
        # eliminate all the pixels with 255 values for all 8 bands - AKA not actually part of the shapefile
        out_image_trimmed = out_image_trimmed[:,~np.all(out_image_trimmed == 255, axis=0)]
        # reshape the array to [pixel count, bands]
        out_image_reshaped = out_image_trimmed.reshape(-1, band_count)
        # append the labels to the y array
        y = np.append(y,[shapefile[Class][index]] * out_image_reshaped.shape[0]) 
        # stack the pizels onto the pixel array
        X = np.vstack((X,out_image_reshaped))
    return X,y
src = full_dataset
X,y = setup_features_and_target(src,geoms,shapefile)
X_test,y_test = setup_features_and_target(src,geoms_test,shapefile_test)
# What are our classification labels?
labels = np.unique(shapefile[Class])
print('The training data include {n} classes: {classes}\n'.format(n=labels.size, 
                                                                classes=labels))
labels_test = np.unique(shapefile_test[Class])
print('The testing data include {n} classes: {classes}\n'.format(n=labels_test.size, 
                                                                classes=labels_test))

# We will need a "X" matrix containing our features, and a "y" array containing our labels
print('Our X matrix is sized: {sz}'.format(sz=X.shape))
print('Our y array is sized: {sz}'.format(sz=y.shape))

print('Our X_test matrix is sized: {sz}'.format(sz=X_test.shape))
print('Our y_test array is sized: {sz}'.format(sz=y_test.shape))

fig, ax = plt.subplots(1,3, figsize=[20,8])

# numbers 1-4
band_count = np.arange(1,5)

classes = np.unique(y)
for class_type in classes:
    band_intensity = np.mean(X[y==class_type, :], axis=0)
    ax[0].plot(band_count, band_intensity, label=class_type)
    ax[1].plot(band_count, band_intensity, label=class_type)
    ax[2].plot(band_count, band_intensity, label=class_type)
# plot them as lines

# Add some axis labels
ax[0].set_xlabel('Band #')
ax[0].set_ylabel('Reflectance Value')
ax[1].set_ylabel('Reflectance Value')
ax[1].set_xlabel('Band #')
ax[2].set_ylabel('Reflectance Value')
ax[2].set_xlabel('Band #')
#ax[0].set_ylim(32,38)
ax[1].set_ylim(32,38)
ax[2].set_ylim(70,140)
#ax.set
ax[1].legend(loc="upper right")
# Add a title
ax[0].set_title('Band Intensities Full Overview')
ax[1].set_title('Band Intensities Lower Ref Subset')
ax[2].set_title('Band Intensities Higher Ref Subset')
def str_class_to_int(class_array):
    class_array[class_array == 'Territórios artificializados'] = 0
    class_array[class_array == 'Agricultura'] = 1
    class_array[class_array == 'Floresta'] = 2
    class_array[class_array == 'Vegetação natural'] = 3
    class_array[class_array == 'Espaços descobertos ou com vegetação esparsa'] = 4
    class_array[class_array == 'Zonas húmidas'] = 5
    class_array[class_array == 'Corpos de água'] = 6
    return(class_array.astype(int))
def str_class_to_int_v2(class_array,class_test):
    uniques = np.unique(class_array)
    i=0
    for class_name in uniques :
        class_array[class_array == class_name ] = i
        class_test[class_test == class_name ] = i
        i+=1
    return(class_array.astype(int) , class_test.astype(int) ) 
print(np.unique(y))
print(np.unique(y_test))
y = str_class_to_int(y)
y_test = str_class_to_int(y_test)
print(np.unique(y))
print(np.unique(y_test))<jupyter_output>['Agricultura' 'Corpos de água'
 'Espaços descobertos ou com vegetação esparsa' 'Floresta'
 'Territórios artificializados' 'Vegetação natural']
['Agricultura' 'Corpos de água']
[0 1 2 3 4 6]
[1 6]<jupyter_text>Supervised Classification Algorithm<jupyter_code>from sklearn.svm import LinearSVC
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsClassifier
#model = LinearSVC(C=1,verbose=1)
model = KNeighborsClassifier()
model.fit(X, y)
score = model.score(X_test, y_test)
score
with rasterio.open(img_fp) as src:
    # may need to reduce this image size if your kernel crashes, takes a lot of memory
    img = src.read()[:, xmin_test : xmax_test ,  ymin_test : ymax_test ]

# Take our full image and reshape into long 2d array (nrow * ncol, nband) for classification
print(img.shape)
reshaped_img = reshape_as_image(img)
print(reshaped_img.shape)
class_prediction = model.predict(reshaped_img.reshape(-1, 4))
before = class_prediction.shape
# Reshape our classification map back into a 2D matrix so we can visualize it
class_prediction = class_prediction.reshape(reshaped_img[:, :, 0].shape)
class_prediction.shape
"""
reshaped_img.reshape(-1, 4).shape
real_img = y_test[:250000] #fixme
print(before)
print(real_img.shape)
print(y_test.shape)
# Reshape our classification map back into a 2D matrix so we can visualize it
real_class = real_img.reshape(reshaped_img[:, :, 0].shape) #fixme
real_class.shape
"""
real_class = real_class.astype(int)
print(real_class)
class_prediction = class_prediction.astype(int)
class_prediction<jupyter_output><empty_output><jupyter_text>Visualizing results<jupyter_code># find the highest pixel value in the prediction image
n = int(np.max(class_prediction))

# next setup a colormap for our map
colors = dict((
    (0, (60, 60, 60, 255)),     # Grey - Territórios artificializados
    (1, (20,230,20, 255)),      # Light Green - Agricultura
    (2, (5, 80, 5, 255)),       # Dark Green - Floresta
    (3, (170, 200, 20, 255)),   # Yellow - Vegetação natural
    (4, (50, 20, 0, 255)),      # Brown - Espaços descobertos ou com vegetação esparsa
    (5, (0, 250, 250, 255)),    # Light Blue - Zonas húmidas
    (6, (0, 0, 200, 255)),      # Dark Blue - Corpos de água
))

# Put 0 - 255 as float 0 - 1
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 
                (255, 255, 255, 0) for key in range(0, n+1)]

cmap = plt.matplotlib.colors.ListedColormap(index_colors, 'Classification', n+1)
def discrete_matshow(data):
    #get discrete colormap
    cmap = plt.get_cmap('RdBu', np.max(data)-np.min(data)+1)
    # set limits .5 outside true range
    mat = plt.matshow(data,cmap=cmap,vmin = np.min(data)-.5, vmax = np.max(data)+.5)
    #tell the colorbar to tick at integers
    cax = plt.colorbar(mat, ticks=np.arange(np.min(data),np.max(data)+1))

fig, axs = plt.subplots(1,2,figsize=(20,20))

img_stretched = color_stretch(reshaped_img, [2, 1, 0])
axs[0].imshow(img_stretched)

prediction_plot= axs[1].imshow(class_prediction, cmap=cmap, interpolation='none')
fig.colorbar(prediction_plot, ax=axs[1],ticks=range(0,6),fraction=0.046, pad=0.04)

#real_plot = axs[2].imshow(real_class, cmap=cmap, interpolation='none')
#fig.colorbar(real_plot, ax=axs[2],ticks=range(0,6),fraction=0.046, pad=0.04)

fig.show()
fig, axs = plt.subplots(1,2,figsize=(15,15))

img_stretched = color_stretch(reshaped_img, [2, 1, 0])

axs[0].imshow(img_stretched[0:100, 0:100])

prediction_plot = axs[1].imshow(class_prediction[0:100, 0:100], cmap=cmap, interpolation='none')
fig.colorbar(prediction_plot, ax=axs[1],ticks=range(0,4),fraction=0.046, pad=0.04)

fig.show()
fig, axs = plt.subplots(2,2,figsize=(15,7))

reshaped_img =reshape_as_image( img )
img_stretched = color_stretch(reshaped_img, [2, 1, 0])
axs[0,0].set_title("True Color")
axs[0,0].imshow(img_stretched)

axs[0,1].set_title("Class Prediction")
predictions_plot = axs[0,1].imshow(class_prediction, cmap=cmap, interpolation='none')
fig.colorbar(predictions_plot, ax=axs[0,1],ticks=range(0,6))

nwdi_plot = axs[1,0].imshow(ndwi, cmap="RdYlBu")
axs[1,0].set_title("NDWI")
fig.colorbar(nwdi_plot, ax=axs[1,0])

ndvi_plot = axs[1,1].imshow(ndvi, cmap="RdYlGn")
axs[1,1].set_title("NDVI")
fig.colorbar(ndvi_plot, ax=axs[1,1])

plt.show()
alert() # audio alert<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13003529900332225, 'fraction_numerical': 0.02616279069767442, 'mean_word_length': 4.898652786282915, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 40, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '233235', 'license': 'MIT', 'n_tokens_mistral': 7540, 'n_tokens_neox': 7067, 'n_words': 1594, 'path': 'Supervised_Classification_of_IMG_DATA_1.ipynb', 'repo_name': 'SimaoBarbosa/Land_cover_classification', 'size': '19264'}