<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Read and plot images2&nbsp;&nbsp;Plot individual band3&nbsp;&nbsp;Read and reproject polygons4&nbsp;&nbsp;Define function to extract statistics for each polygon5&nbsp;&nbsp;Iterate all polygons and derive data into a panda dataframe<jupyter_code>import os
import rasterio as rio
from rasterio.plot import show
import pandas as pd
import geopandas as gpd
import glob
%matplotlib inline
from matplotlib import mlab 
import matplotlib.pyplot as plt
import numpy as np
from rasterio import mask
import fiona
import json
import multiprocessing as mp
from time import time
import warnings
from numpy import inf
warnings.filterwarnings("ignore")
# os.getcwd()<jupyter_output><empty_output><jupyter_text>Read and plot images<jupyter_code>tifs = glob.glob('*.tif')
#tifs
img = rio.open('**.tif')
show(img)<jupyter_output><empty_output><jupyter_text>Plot individual band<jupyter_code>arys = img.read()
arys = arys.astype('uint8')
arys2 = np.moveaxis(arys, 0, -1)
print(arys2.shape)
b1 = arys2[:,:,0]
b2 = arys2[:,:,1]
b3 = arys2[:,:,2]
b4 = arys2[:,:,3]
# band 1 - red
plt.imshow(b1)
# band 2 - green
plt.imshow(b2)
# band 3 - blue
plt.imshow(b3)
# band 4 - empty
plt.imshow(b4)<jupyter_output><empty_output><jupyter_text>Read and reproject polygons<jupyter_code>ply = gpd.read_file('*.shp')
print(ply.shape)
ply.plot()
print(ply.crs)
ply2 = ply.to_crs('EPSG:32616')
#ply2.to_file("polys.shp")<jupyter_output>{'init': 'epsg:4326'}<jupyter_text>Define function to extract statistics for each polygon<jupyter_code>def extract_plot(i, image,  ply):
    coords = [json.loads(ply.iloc[[i]].to_json())["features"][0]["geometry"]]
    out_image, out_transform  = mask.mask(image, coords, crop=True)
    arys2 = np.moveaxis(out_image, 0, -1)
    
    r = arys2[:,:,0]
    g = arys2[:,:,1]
    b = arys2[:,:,2]
    r= r.astype(np.int32)
    g= g.astype(np.int32)
    b= b.astype(np.int32)
    #nir = arys2[:,:,3]

    #nv = (nir-r)/(nir+r)
    bi = np.sqrt((r*r + g*g +b*b)/3)
    sci = (r-g)/(r+g)
    gli = (2*g-r-b)/(2*g+r+b)
    hi = (2*r-g-b)/(g-b)

    ngrdi = (g-r)/(g+r)
    si = (r-b)/(r+b)
    vari = (g-r)/(g+r-b)
    bgi = b/g
    hue = np.arctan(2*(b-g-r)/(30.5*(g-r)))
    
    vis = np.dstack((bi, sci, gli, hi, ngrdi, si, vari, hue, bgi))
    vis[np.isinf(vis)] = np.nan
    vis = vis.astype('float64')

    
    mds = np.nanmedian(vis, axis=(0,1))
    stds = np.nanstd(vis, axis=(0,1))
    means = np.nanmean(vis, axis=(0,1))

    vi_nms = ['bi', 'sci', 'gli', 'hi', 'ngrdi', 'si', 'vari', 'hue', 'bgi']
    vi_nms_md = ['median_' + nm for nm in vi_nms]
    vi_nms_mn = ['mean_' + nm for nm in vi_nms]
    vi_nms_std = ['std_' + nm for nm in vi_nms]

    mds_df = pd.DataFrame(mds).T
    stds_df = pd.DataFrame(stds).T
    means_df = pd.DataFrame(means).T

    mds_df.columns = vi_nms_md
    stds_df.columns = vi_nms_std
    means_df.columns = vi_nms_mn

    all_df = pd.concat([means_df, mds_df, stds_df],axis=1)
    
    return all_df<jupyter_output><empty_output><jupyter_text>Iterate all polygons and derive data into a panda dataframe<jupyter_code>import time
t0 = time.time()

results_df = pd.DataFrame()

for i in range(100):
    vi_values = extract_plot( i, img, ply2)
    results_df = pd.concat([results_df, vi_values], axis=0)
    
    print('Progress is: ' + str(i*100/100) + "%")

print("Time consumed is: " + str(time.time() - t0))
results_df<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16647298082510167, 'fraction_numerical': 0.025857059848925043, 'mean_word_length': 4.1852409638554215, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '441892', 'license': 'MIT', 'n_tokens_mistral': 1535, 'n_tokens_neox': 1471, 'n_words': 302, 'path': 'code/Extract raster value by polygons.ipynb', 'repo_name': 'wanwanliang/Spatial-Aanalysis-with-Python', 'size': '3442'}