<jupyter_start><jupyter_text>GGS416 Satellite Image AnalysisIn this tutorial we are going to explore `numpy` and `rasterio`.As the `numpy` library provides support for large, multi-dimensional arrays and matrices, the package is widely used by a variety (in fact, most) other python packages, including: `pandas`, `scipy`, `matplotlib`, `scikit-learn`, `scikit-image` etc.See here for more details: https://en.wikipedia.org/wiki/NumPyLet's begin by loading an image we have already downloaded:<jupyter_code># First we need to import rasterio (Remember: it stands for raster input output - rasterio)
import rasterio

# Then we can open our image file and create a rasterio object
my_image = rasterio.open("../week3/20190321_174348_0f1a_3B_AnalyticMS.tif")
my_image<jupyter_output><empty_output><jupyter_text>As this is a 4-band PlanetScope image, we are then able to parse out these image bands to new variables.<jupyter_code># So we can use the rasterio read function and allocate our bands to blue, green, red and near-infrared
blue, green, red, nir = my_image.read()

# Each one is a numpy array of our data:
blue<jupyter_output><empty_output><jupyter_text>Each `numpy` array is essentially a grid of different pixels with values.We can confirm this is a `numpy` array by checking the data type!<jupyter_code># Bands are stored as Numpy arrays.
print(type(blue))<jupyter_output><class 'numpy.ndarray'><jupyter_text>"A raster graphic is a mechanism that represents a two-dimensional image as a rectangular matrix or grid of square pixels. A raster is technically characterized by the width and height of the image in pixels and by the number of bits per pixel."See here for more info on raster structures: https://en.wikipedia.org/wiki/Raster_graphicsSo let's look at this Digital Elevation Model raster below for Humboldt County in Northern California:<jupyter_code>Image("images/RasterPixels.png")

# Credit to 
# https://gsp.humboldt.edu/olm/Lessons/GIS/08%20Rasters/RasterDataModels3.html<jupyter_output><empty_output><jupyter_text>This should give you an idea of how these raster layers are composed. Now let's take a look at satellite image raster pixels:<jupyter_code>Image("images/pixels.png")

# Credit to 
# https://mapschool.io/<jupyter_output><empty_output><jupyter_text>The difference with color satellite imagery is that our raster pixels are indicating color shades.Below is a stylized example of how pixels showing comparable colors have similar underlying data:<jupyter_code>Image("images/Pixels2.png")

# Credit to 
# https://github.com/sarasafavi/remote-sensing-with-python/tree/master/notebooks<jupyter_output><empty_output><jupyter_text>Thus, a single raster band has two dimensions consisting of rows and columns.<jupyter_code># We can confirm this using the ndim function.
# ndim stands for number of dimensions
print(blue.ndim)<jupyter_output>2<jupyter_text>We can also explore the type of underlying data that we have stored in these pixel grids:<jupyter_code># So here we have 16-bit unsigned integers
# These are integers which can range between 0 and 65,535
print(blue.dtype)<jupyter_output>uint16<jupyter_text>For all the pixels in our data, we are able to use the `min()` and `max()` functions to find the range of values.- `min()` represents the minimum data value in a given array. - `max()` represents the maximum data value in a given array.<jupyter_code># Let's print the minimum and maximum values
list_of_our_bands = [blue, green, red, nir]
for band in list_of_our_bands:
    print("min {min} max {max}".format(min=band.min(), max=band.max()))<jupyter_output>min 0 max 30256
min 0 max 20671
min 0 max 27564
min 0 max 33932<jupyter_text>We can extract information from the raster layer using indexing.For example, here we can try to obtain the values for a pixel which is 5 km east and 5 km south of the top left corner. Let's get the x and y coordinates first. Remember our Planet image is in the Universal Transverse Mercator (UTM) projection, and here we are measuring in meters (so not in degrees or feet, for example).<jupyter_code>import math

# World coordinates for the desired pixel.
x_coord = my_image.bounds.left + 5000
y_coord = my_image.bounds.top - 5000
print('Our x and y cordinates are: {}, {}'.format(x_coord, y_coord))<jupyter_output>Our x and y cordinates are: 549491.0, 4186009.0<jupyter_text>Next, we can use the `index()` function for our `rasterio` object. This converts for our world coordinates in UTM into row and column indexes.This may not be a precise transform, so the `op` argument allows us to control for how values are rounded.<jupyter_code># So we give the index function our coordinates
# math.floor round numbers down to the nearest integer
col, row = my_image.index(x_coord, y_coord, op=math.floor)

# Let's make sure we have integers
col = int(col)
row = int(row)

print(col, row)<jupyter_output>1666 1666<jupyter_text>Now we have our column and row indexes, we can use them to index into our `numpy` array raster data.This demonstrates that we have an integer number of 2344 for red, 3658 for green, 4445 for blue, and 969 for NIR.<jupyter_code># Now let's look at the value of each band at this pixel
print("Red: {}".format(red[row, col]))
print("Green: {}".format(green[row, col]))
print("Blue: {}".format(blue[row, col]))
print("NIR: {}".format(nir[row, col]))<jupyter_output>Red: 2344
Green: 3658
Blue: 4445
NIR: 969
{'fraction_non_alphanumeric': 0.08486196034834168, 'fraction_numerical': 0.023346303501945526, 'mean_word_length': 5.519323671497585, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 43, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '531930', 'license': 'MIT', 'n_tokens_mistral': 1768, 'n_tokens_neox': 1635, 'n_words': 682, 'path': 'notebooks/week4/4.02_intro_to_numpy_and_rasterio.ipynb', 'repo_name': 'edwardoughton/Satellite-Image-Analysis', 'size': '5397'}