<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Inspecting Satellite Imagery using Rasterio
A first look at satellite data with Python
At this point, you will have learned different ways of searching for, filtering, and downloading satellite imagery. Now let's use one of these acquired datasets and dig into it a bit with Python.
Here we're going to use a Python library called Rasterio
Step1:  Basic details
What can we learn about this satellite image using just Python?
Step2:  Bands
So far, we haven't done too much geospatial-raster-specific work yet. Since we know we're inspecting a multispectral satellite image, let's see what we can learn about its bands.
Step3:  Because we know we're look at a PlanetScope 4-band analytic satellite image, we can define the bands by their order
Step4:  Pixels
In a raster dataset, each pixel has a value. Pixels are arranged in a grid, and pixels representing equivalent data have the same value
<ASSISTANT_TASK:>
Python Code:
from __future__ import division
import math
import rasterio
# This notebook explores a single 4 band (blue, green, red, NIR) PlanetScope scene in a UTM projection.
image_file = "example.tif"
satdat = rasterio.open(image_file)
Explanation: Inspecting Satellite Imagery using Rasterio
A first look at satellite data with Python
At this point, you will have learned different ways of searching for, filtering, and downloading satellite imagery. Now let's use one of these acquired datasets and dig into it a bit with Python.
Here we're going to use a Python library called Rasterio: you may be familiar with it already, or perhaps with the related C library, GDAL. If you've used Numpy before, working with Rasterio will feel very familiar.
End of explanation
# Minimum bounding box in projected units
print(satdat.bounds)
# Get dimensions, in map units (using the example GeoTIFF, that's meters)
width_in_projected_units = satdat.bounds.right - satdat.bounds.left
height_in_projected_units = satdat.bounds.top - satdat.bounds.bottom
print("Width: {}, Height: {}".format(width_in_projected_units, height_in_projected_units))
# Number of rows and columns.
print("Rows: {}, Columns: {}".format(satdat.height, satdat.width))
# This dataset's projection uses meters as distance units.  What are the dimensions of a single pixel in meters?
xres = (satdat.bounds.right - satdat.bounds.left) / satdat.width
yres = (satdat.bounds.top - satdat.bounds.bottom) / satdat.height
print(xres, yres)
print("Are the pixels square: {}".format(xres == yres))
# Get coordinate reference system
satdat.crs
# Convert pixel coordinates to world coordinates.
# Upper left pixel
row_min = 0
col_min = 0
# Lower right pixel.  Rows and columns are zero indexing.
row_max = satdat.height - 1
col_max = satdat.width - 1
# Transform coordinates with the dataset's affine transformation.
topleft = satdat.transform * (row_min, col_min)
botright = satdat.transform * (row_max, col_max)
print("Top left corner coordinates: {}".format(topleft))
print("Bottom right corner coordinates: {}".format(botright))
# All of the metadata required to create an image of the same dimensions, datatype, format, etc. is stored in
# one location.
print(satdat.meta)
Explanation: Basic details
What can we learn about this satellite image using just Python?
End of explanation
# The dataset reports a band count.
print(satdat.count)
# And provides a sequence of band indexes.  These are one indexing, not zero indexing like Numpy arrays.
print(satdat.indexes)
Explanation: Bands
So far, we haven't done too much geospatial-raster-specific work yet. Since we know we're inspecting a multispectral satellite image, let's see what we can learn about its bands.
End of explanation
# PlanetScope 4-band band order: BGRN
blue, green, red, nir = satdat.read()
# Or the slightly less efficient:
#    blue = satdat.read(1)
#    green = satdat.read(2)
#    red = satdat.read(3)
#    nir = satdat.read(4)
# Or read the entire dataset into a single 3D array:
#    data = satdat.read()
Explanation: Because we know we're look at a PlanetScope 4-band analytic satellite image, we can define the bands by their order:
End of explanation
# Bands are stored as Numpy arrays.
print(type(blue))
# How many dimensions would a single raster band have?  Two dimensions: rows and columns.
print(blue.ndim)
# Glimpse at the band's values and datatype.
print(blue)
print(blue.dtype)
# Output a min & max pixel value in each band.
for bidx in satdat.indexes:
    data = satdat.read(bidx)
    print("Band {bidx} min {min} max {max}".format(bidx=bidx, min=data.min(), max=data.max()))
# And an overall min/max for the entire dataset.
    
data = satdat.read()
print("Overall min/max: {} {}".format(data.min(), data.max()))
# Let's grab the pixel 2km east and 2km south of the upper left corner
# World coordinates for the desired pixel.
x_coord = satdat.bounds.left - 2000
y_coord = satdat.bounds.top + 2000
# Convert world coordinates to pixel.  World coordinates may not transform precisely to row and column indexes,
# but a Numpy array can only be indexed by integer values.  The 'op' parameter for 'satdat.index()' determines
# how the transformed values are rounded.  In some cases any point falling within a pixel should be considered
# contained, and in other cases only points falling within one portion of the pixels hould be considered contained.
# The 'op' parameter lets users make this decision on their own.  The values must still be cast to integers.
col, row = satdat.index(x_coord, y_coord, op=math.floor)
col = int(col)
row = int(row)
# 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]))
Explanation: Pixels
In a raster dataset, each pixel has a value. Pixels are arranged in a grid, and pixels representing equivalent data have the same value:
End of explanation
<END_TASK>