<jupyter_start><jupyter_text><jupyter_code>#https://www.earthdatascience.org/courses/use-data-open-source-python/multispectral-remote-sensing/vegetation-indices-in-python/calculate-NDVI-python/
#Sentinel 2 Use Handbook. https://sentinels.copernicus.eu/documents/247904/685211/Sentinel-2_User_Handbook


!pip install rioxarray

!pip install geopandas
import os

import matplotlib.pyplot as plt
import numpy as np
import rioxarray as rxr
import geopandas as gpd


path = '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B01.tif'
#os.chdir(path)
#data_path = os.path.join("/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B01.tif")
data = rxr.open_rasterio(path)
data.shape
!pip install earthpy
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
ep.plot_bands(data,
              title="Bigearthnet Band 1 Raster")
plt.show()
#https://earthpy.readthedocs.io/en/latest/gallery_vignettes/plot_bands_functionality.html
#Stack all bands of BigEarthNet Data sample one band tiff images
import glob
files = glob.glob(os.path.join('/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B*.tif'))
files.sort()
print("Number of Bands",len(files))
print(files)
print(files[0])   # Band1
print(files[1])   # Band2
print(files[10])  # Band12
#array_stack, meta_data = es.stack(path, nodata=-9999)<jupyter_output>Number of Bands 12
['/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B01.tif', '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B02.tif', '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B03.tif', '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B04.tif', '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B05.tif', '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B06.tif', '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B07.tif', '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B08.tif', '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B09.tif', '/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2[...]<jupyter_text>New Section<jupyter_code>print(files[0])
band1= rxr.open_rasterio(files[0])
ep.plot_bands(band1,
              title="Bigearthnet Band 1 Raster")
plt.show()
print("The CRS of this data is:", band1.rio.crs)
#Converting EPSG to Proj4 in Python
# Convert to project string using earthpy
proj4 = et.epsg['32634']
print(proj4)
#Spatial Extent
#You can access the spatial extent using the .bounds() attribute in rasterio.
print(band1.rio.bounds())
#Raster Resolution: area covered by 1 pixel on ground e.g 60m * 60m
# What is the x and y resolution for your raster data?
print(band1.rio.resolution())
print("The nodatavalue of your data is:", band1.rio.nodata)
# How many bands / layers does the object have?
print("Number of bands", band1.rio.count)
print("The shape of your data is:", band1.shape)
print('min value:', np.nanmin(band1))
print('max value:', np.nanmax(band1))
import matplotlib.pyplot as plt
f, ax = plt.subplots()
band1.plot.hist(color="purple")
ax.set(title="Distribution of Raster Cell Values  Band 1 Data",
       xlabel="",
       ylabel="Number of Pixels")
plt.show()
print(files[1])
band2= rxr.open_rasterio(files[1])
ep.plot_bands(band2,
              title="Bigearthnet Band 2 Raster")
plt.show()
print("The CRS of this data is:", band2.rio.crs)
#Converting EPSG to Proj4 in Python
# Convert to project string using earthpy
proj4 = et.epsg['32634']
print(proj4)
#Spatial Extent
#You can access the spatial extent using the .bounds() attribute in rasterio.
print(band2.rio.bounds())
#Raster Resolution: area covered by 1 pixel on ground e.g 60m * 60m
# What is the x and y resolution for your raster data?
print(band2.rio.resolution())
print("The nodatavalue of your data is:", band2.rio.nodata)
# How many bands / layers does the object have?
print("Number of bands", band2.rio.count)
print("The shape of your data is:", band2.shape)
print('min value:', np.nanmin(band2))
print('max value:', np.nanmax(band2))
import matplotlib.pyplot as plt
f, ax = plt.subplots()
band1.plot.hist(color="purple")
ax.set(title="Distribution of Raster Cell Values  Band 2 Data",
       xlabel="",
       ylabel="Number of Pixels")
plt.show()
#https://rasterio.readthedocs.io/en/latest/api/rasterio.plot.html
#rasterio.plot.reshape_as_image(arr)
#Returns the source array reshaped into the order expected by image processing and visualization software (matplotlib, scikit-image, etc) by swapping the axes order from (bands, rows, columns) to (rows, columns, bands)

print('min value:', np.nanmin(data))
print('max value:', np.nanmax(data))
#https://www.earthdatascience.org/courses/use-data-open-source-python/intro-raster-data-python/raster-data-processing/classify-plot-raster-data-in-python/
import matplotlib.pyplot as plt
f, ax = plt.subplots()
data.plot.hist(color="purple")
ax.set(title="Distribution of Raster Cell Values  Data",
       xlabel="",
       ylabel="Number of Pixels")
plt.show()
bins=[0, 100, 200, 250, 275, 300,350]
f, ax = plt.subplots()
data.plot.hist(color="purple",bins=[0, 100, 200, 250, 275, 300,350])
ax.set(title="Distribution of Raster Cell Values  Data",
       xlabel="",
       ylabel="Number of Pixels")
plt.show()
class_bins = [-np.inf,250,275,300,350,+np.inf]

import xarray as xr
data_class = xr.apply_ufunc(np.digitize,
                                     data,
                                     class_bins)
print(data_class.shape)
#data_class = np.array(data_class[0])
import matplotlib.pyplot as plt
f, ax = plt.subplots()
data_class.plot.hist(color="purple")
ax.set(title="Distribution of Raster Cell Values  Data",
       xlabel="",
       ylabel="Number of Pixels")
plt.show()
#https://www.spatialreference.org/ref/epsg/32634/
#/*EPSG:32634

#WGS 84 / UTM zone 34N (Google it)

#WGS84 Bounds: 18.0000, 0.0000, 24.0000, 84.0000
#Projected Bounds: 166021.4431, 0.0000, 833978.5569, 9329005.1825
#Scope: Large and medium scale topographic mapping and engineering survey.
#Last Revised: June 2, 1995
#Area: World - N hemisphere - 18°E to 24°E - by country*/

#Proj4js.defs["EPSG:32634"] = "+proj=utm +zone=34 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
print(files[10])
band12= rxr.open_rasterio(files[10])
ep.plot_bands(band12,
              title="Bigearthnet Band 12 Raster")
plt.show()
print("The CRS of this data is:", band12.rio.crs)
#Converting EPSG to Proj4 in Python
# Convert to project string using earthpy
proj4 = et.epsg['32634']
print(proj4)
#Spatial Extent
#You can access the spatial extent using the .bounds() attribute in rasterio.
print(band12.rio.bounds())
#Raster Resolution: area covered by 1 pixel on ground e.g 60m * 60m
# What is the x and y resolution for your raster data?
print(band12.rio.resolution())
print("The nodatavalue of your data is:", band12.rio.nodata)
# How many bands / layers does the object have?
print("Number of bands", band12.rio.count)
print("The shape of your data is:", band12.shape)
print('min value:', np.nanmin(band12))
print('max value:', np.nanmax(band12))
import matplotlib.pyplot as plt
f, ax = plt.subplots()
band1.plot.hist(color="purple")
ax.set(title="Distribution of Raster Cell Values  Band 12 Data",
       xlabel="",
       ylabel="Number of Pixels")
plt.show()<jupyter_output>/content/drive/MyDrive/Big/S2A_MSIL2A_20170613T101031_0_55/S2A_MSIL2A_20170613T101031_0_55_B12.tif<jupyter_text>New Section : Making Dataframe for min/max values of each bands of 1,2 and 12<jupyter_code>import pandas as pd
df = pd.DataFrame(columns= ['filename','min','max'])
df.head()
import glob
import os
files_batch=[]  #batch of same bands 
min=[]
max=[]
mean_min =[]
mean_max =[]
path = '/content/drive/MyDrive/Big'
os.chdir(path)
dirs = os.listdir()
dirs.sort()
print(dirs)
print(len(dirs))
#remove last element of list
del dirs[0]
print(dirs)
print(len(dirs))
step_size = len(dirs)
# Add batch of band1 tif files to files list
path = '/content/drive/MyDrive/Big'

for i in dirs:
  s = ""
  s = s + path + '/' + str(i) +  '/' +'*01.tif'
  print(s)
  temp = (glob.glob(os.path.join(s)))
  files_batch.append(temp[0])

# Fetch Filenames of band 1

print(files_batch,files_batch[0],len(files_batch)) #Batch of Band 1 files
# Add min/max values of band 1 to min/max list

for i in range(0,step_size):
  band1= rxr.open_rasterio(files_batch[i])
  min.append(np.nanmin(band1))
  max.append(np.nanmax(band1))
print(min)
print(max)
mean_min.append(np.mean(min))
mean_max.append(np.mean(max))
#df['B1_min'] = min
#df['B1_max'] = max
#print(df)
# Add batch of band2 tif files to files list
path = '/content/drive/MyDrive/Big'

for i in dirs:
  s = ""
  s = s + path + '/' + str(i) +  '/' +'*02.tif'
  print(s)
  temp = (glob.glob(os.path.join(s)))
  files_batch.append(temp[0])
print(files_batch)
print(files_batch[len(files_batch)-1], len(files_batch))

# Add min/max values of band 2 to min/max list

for i in range(step_size,2*step_size):
  band2= rxr.open_rasterio(files_batch[i])
  min.append(np.nanmin(band2))
  max.append(np.nanmax(band2))
print(min)
print(max)
mean_min.append(np.mean(min))
mean_max.append(np.mean(max))
# Add batch of band 12 tif files to files list
path = '/content/drive/MyDrive/Big'

for i in dirs:
  s = ""
  s = s + path + '/' + str(i) +  '/' +'*12.tif'
  print(s)
  temp = (glob.glob(os.path.join(s)))
  files_batch.append(temp[0])
print(files_batch)
print(files_batch[len(files_batch)-1], len(files_batch))

# Add min/max values of band 12 to min/max list

for i in range(2*step_size,3*step_size):
  band2= rxr.open_rasterio(files_batch[i])
  min.append(np.nanmin(band2))
  max.append(np.nanmax(band2))
print(min)
print(max)
mean_min.append(np.mean(min))
mean_max.append(np.mean(max))<jupyter_output>[267, 272, 294, 280, 245, 276, 123, 156, 198, 164, 139, 202, 92, 91, 95, 104, 119, 115]
[456, 308, 655, 327, 302, 328, 1215, 367, 985, 567, 515, 456, 149, 144, 1606, 153, 191, 170]<jupyter_text>Add files and min/max lists to dataframe<jupyter_code>print(files_batch)
df['filename'] = files_batch
df['min'] = min
df['max'] = max
df.head()
#print means of min and max values for each band 1 2 and 12
print(mean_min)
print(mean_max)
# Plot histogram 
import matplotlib.pyplot as plt
x=np.array(min)
y=np.array(max)
plt.bar(x,y,align='center') # A bar chart
plt.xlabel('Min')
plt.ylabel('Max')
plt.show()
# Plot histogram for mean min and mean max
import matplotlib.pyplot as plt
x=np.array(mean_min)
y=np.array(mean_max)
plt.bar(x,y,align='center') # A bar chart
plt.xlabel('Mean_Min')
plt.ylabel('Mean_Max')
plt.show()<jupyter_output><empty_output><jupyter_text>**USE RASTERIO module to open Raster images and read it to Array**<jupyter_code>band1 = np.array(band1)
band1.shape
print(files)
band2= rxr.open_rasterio(files[1])
band2 = np.array(band2)
band2.shape
band12 = np.array(band12)
band12.shape
print(df['filename'])
files_bands = []
files_bands = df['filename']
print(files_bands[0:6])
# Reading raster geotif files 
#https://automating-gis-processes.github.io/CSC18/lessons/L6/reading-raster.html
import rasterio 
band1_batch = files_bands[0:6]
print(band1_batch[0])
band1_raster = rasterio.open(band1_batch[0])
print(type(band1_raster))
#Projection
print(band1_raster.crs)
#Affine transform (how raster is scaled, rotated, skewed, and/or translated
band1_raster.transform
band1_raster.meta
#reading raster to array

band1_array = band1_raster.read()
print(band1_array)
stats = []

for band in band1_array:
  stats.append({
      'mean' :  band.mean(),
      'min'  :  band.min(),
      'max'  :  band.max(),
      'median': np.median(band)
  })
print(stats)<jupyter_output>[{'mean': 284.155, 'min': 267, 'max': 456, 'median': 282.0}]<jupyter_text>Read all Band1 files and find mean of all 6 Forest class Band1 data<jupyter_code>print(df['filename'])
files_bands = []
files_bands = df['filename']
print(files_bands[0:6])
# Reading raster geotif files using Rasterio
#https://automating-gis-processes.github.io/CSC18/lessons/L6/reading-raster.html
import rasterio 
band1_batch = files_bands[0:6]
print(band1_batch[0])
band1_array=[]
for i in band1_batch:
  band1_raster = rasterio.open(i)
  band1_array.append(band1_raster.read())
band1_mean=[]
band1_min = []
band1_max = []
print(len(band1_array))
for i in band1_array:
    for band in i:
        band1_mean.append(band.mean())
        band1_min.append(band.min())
        band1_max.append(band.max())
print("Band 1 stat for 6 images is :------>")
print(band1_mean)
print(band1_min)
print(band1_max)
# Stat for band 2 images 
band2_batch = files_bands[6:12]
print(band2_batch)
band2_array=[]
for i in band2_batch:
  band2_raster = rasterio.open(i)
  band2_array.append(band2_raster.read())
band2_mean=[]
band2_min = []
band2_max = []
print(len(band2_array))
for i in band2_array:
    for band in i:
        band2_mean.append(band.mean())
        band2_min.append(band.min())
        band2_max.append(band.max())
print("Band 2 stat for 6 images is :------>")
print(band2_mean)
print(band2_min)
print(band2_max)
# Stat for band 12 images 
band12_batch = files_bands[12:18]
print(band12_batch)
band12_array=[]
for i in band12_batch:
  band12_raster = rasterio.open(i)
  band12_array.append(band12_raster.read())
band12_mean=[]
band12_min = []
band12_max = []
print(len(band12_array))
for i in band12_array:
    for band in i:
        band12_mean.append(band.mean())
        band12_min.append(band.min())
        band12_max.append(band.max())
print("Band 12 stat for 6 images is :------>")
print(band12_mean)
print(band12_min)
print(band12_max)
y=np.array(band1_mean)
x=(1,2,3,4,5,6)
plt.bar(x,y,align='center') 
plt.axis([0, 6, 100, 600])
plt.xlabel('Bands')
plt.ylabel('Commulative Mean')
plt.show()
y=np.array(band2_mean)
x=(1,2,3,4,5,6)
plt.bar(x,y,align='center') 
plt.axis([0, 6, 100, 600])
plt.xlabel('Bands')
plt.ylabel('Commulative Mean')
plt.show()
y=np.array(band12_mean)
x=(1,2,3,4,5,6)
plt.bar(x,y,align='center') 
plt.axis([0, 6, 100, 600])
plt.xlabel('Bands')
plt.ylabel('Commulative Mean')
plt.show()
df = pd.DataFrame(columns = ['mean_band1','mean_band2','mean_band12'])
df['mean_band1'] = np.array(band1_mean)
df['mean_band2'] = np.array(band2_mean)
df['mean_band12'] = np.array(band12_mean)
df
df.plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1427494000685636, 'fraction_numerical': 0.08049365786767226, 'mean_word_length': 5.238665526090676, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 24, 'https://': 8, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '14191', 'license': 'MIT', 'n_tokens_mistral': 6526, 'n_tokens_neox': 5740, 'n_words': 1225, 'path': 'Understand_band_data_info_using_histogram_and_classifying_pixel_values.ipynb', 'repo_name': 'BRIJNANDA1979/CNN-Sentinel', 'size': '14585'}