Raster Fundamentals Practice#

Part 1#

For the files in this practice we will be using the raster 'subset_f180628t01p00r02_corr_v1k1_img' in the data folder.

Opening a File#

  1. Write a statement to open up the data file 'subset_f180628t01p00r02_corr_v1k1_img', which is located in the data folder in the same folder as this file. Print out the data in the first band of the raster.

  1. Use pyplot to get a quick visual of band 100 of the data file.

Inspecting Metadata#

  1. What is the crs of the raster? (It’s alright if you don’t know what it means yet)

  1. How many rows does this raster have? How many columns? How many bands?

  1. What is the file format?

Working with data#

  1. Read in all of the bands of the raster and save the data to a variable called reflectance. What is the data structure of the reflectance variable?

  1. Get the value of the data at band 110, row 200, column 27

  1. Get all the data for band 20

  1. One common need in remote sensing is to take a z transect. This means that you are taking all from a single pixel location across all of the bands. Pick any location on the raster and write out the index needed to get a z transect for that location. (Bonus question: what does the index look like for an x profile (across a column)? What about a y profile (across a row)?)

  1. What is the mean value of band 23?

Saving a file#

  1. Take one band of data and save it to its own file. Be sure to update the metadata to set count equal to 1 when you save the file.

Part 2#

import numpy as np

Question 1#

The wild looking names of many files is usually a shorthand that includes some metadata about the raster. Use the AVIRIS documentation to figure out the following pieces of metadata from the name of our raster, f180628t01p00r02_corr_v1k1_img.

  • day, month and year of the flight

  • flight run number

Question 2#

Open the raster file. Read all the bands to a variable called reflectance. Assign the src.meta object to a variable called metadata. Use the the metadata dictionary and the shape of the reflectance numpy array to write 3 boolean statements confirming that the count, width, and height in the metadata accurately line up with the .shape of the reflectance array.

Question 3#

np.where() is a really useful function, but the use patterns aren’t always intuitive. In this problem we are going to walk through using np.where().

In this question we will use the following example array:

example_array = np.array([20, 18, 16, 14, 12, 10, 8, 6, 4, 2, 0])
print(example_array)

[20 18 16 14 12 10  8  6  4  2  0]

A) One way to use np.where is to give it a condition. In the following example we are looking for where in our example_array there are values less than 8.

np.where(example_array < 8)

(array([ 7,  8,  9, 10]),)

The result of the array is not the actual values which are less than 8, but the indexes of the values less than 8. So we can check and see that each of the values at those indexes is indeed less than 8:

example_array[7]  # Change 7 to 8, 9, or 10 to see the other values less than 8

6

Question A) Write a condidtional in a np.where to find the values that are greater than 11

B) Since we get as a result all the indexes of the array, if we want the values out of the original array we can index the original array with our np.where() statement, like this:

less_than_8_indexes = np.where(example_array < 8)
example_array[less_than_8_indexes]

array([6, 4, 2, 0])

Question B) Write a statement which returns the values of example_array that are less than or equal to 13.

C) In multi-dimensional arrays the concept stays the same, but we get two sets of indices back from the np.where() statement, because we now have two axes in our array.

example_array_2d = np.array([[2, 4, 6], [8, 10, 12], [14, 16, 18], [20, 22, 24]])
print(example_array_2d)

[[ 2  4  6]
 [ 8 10 12]
 [14 16 18]
 [20 22 24]]

np.where(example_array_2d < 16)

(array([0, 0, 0, 1, 1, 1, 2]), array([0, 1, 2, 0, 1, 2, 0]))

The first array has the indexes of the rows and the second array has the indexes of the columns. Together they make the indexes.

Indexes of numbers below 16, read in order from the output of the np.where() statement:

  • [0,0] (value of 2)

  • [0,1] (value of 4)

  • [0,2] (value of 6)

  • [1,0] (value of 8)

  • [1,1] (value of 10)

  • [1,2] (value of 12)

  • [2,0] (value of 14)

If you’re a kinestetic learner like I am, try these indexes out in the cell below to confirm.

# Change the index value to each of the index combinations from the `np.where` statement
example_array_2d[1,2]

12

Question C) Use np.where to find the 4 index locations of values in example_array_2d that have values greater than 17. Go back and look at example_array_2d to make sure those locations make sense.

Question 4#

One calculation often done in remote sensing is NDVI - the Normalized Difference Vegetation Index. The formula for that is:

NDVI = (r857 - r645) / (r857 + r645)

Where r857 refers to the band with wavelength 857 (call it NIR, or near infrared) and r645 refers to the band of wavelength 645 (red band).

A) The following cell contains an numpy array that I copied from the header file that specifies the wavelengths of each of the bands of the AVIRIS file. The closest bands to the wavelengths we need are 859.64709 and 647.97357. Use np.where() to the index number of each of those wavelengths in the wavelengths array.

wavelengths = np.array([365.92981, 375.59399, 385.26254, 394.93552, 404.61288, 414.29462, 
                        423.98077, 433.6713, 443.36621, 453.06552, 462.7692, 472.47729, 
                        482.18976, 491.90665, 501.6279, 511.35355, 521.08356, 530.81799, 
                        540.55682, 550.30005, 560.04767, 569.79962, 579.55603, 589.31677, 
                        599.08191, 608.8515, 618.62543, 628.40375, 638.18646, 647.97357, 657.76508, 
                        667.56097, 654.7923, 664.59937, 674.40125, 684.19794, 693.98938, 703.77563, 
                        713.55664 , 723.33252 , 733.10309 , 742.86853 , 752.62872 , 762.38373 , 772.13348 , 781.87805 , 791.61743 , 801.35156 , 811.08051 , 820.80426 , 830.52277 , 840.23608 , 849.94415 , 859.64709 , 869.34479 , 879.03723 , 888.72449 , 898.40656 , 908.08337 , 917.75507 , 927.42145 , 937.0827 , 946.73871 , 956.38947 , 966.0351 , 975.67548 , 985.31061 , 994.94055 , 1004.5653 , 1014.18488 , 1023.79919 , 1033.40833 , 1043.01221 , 1052.61096 , 1062.20447 , 1071.79272 , 1081.37573 , 1090.95361 , 1100.52625 , 1110.09375 , 1119.65601 , 1129.21301 , 1138.76477 , 1148.3114 , 1157.85278 , 1167.38904 , 1176.91992 , 1186.4458 , 1195.96631 , 1205.48169 , 1214.99182 , 1224.4967 , 1233.99646 , 1243.49097 , 1252.98022 , 1262.46436 , 1252.77295 , 1262.74573 , 1272.71826 , 1282.69055 , 1292.66248 , 1302.63428 , 1312.60583 , 1322.57715 , 1332.5481 , 1342.51892 , 1352.4895 , 1362.45984 , 1372.42993 , 1382.39978 , 1392.36938 , 1402.33875 , 1412.30786 , 1422.27673 , 1432.24536 , 1442.21375 , 1452.18188 , 1462.1499 , 1472.11755 , 1482.08496 , 1492.05212 , 1502.01904 , 1511.98584 , 1521.95227 , 1531.91846 , 1541.88452 , 1551.85022 , 1561.81567 , 1571.78101 , 1581.74597 , 1591.71082 , 1601.67529 , 1611.63965 , 1621.60364 , 1631.5675 , 1641.53101 , 1651.49438 , 1661.45752 , 1671.42029 , 1681.38293 , 1691.34534 , 1701.30737 , 1711.26929 , 1721.23096 , 1731.19238 , 1741.15344 , 1751.11438 , 1761.07507 , 1771.03552 , 1780.99573 , 1790.95569 , 1800.91541 , 1810.87488 , 1820.83411 , 1830.79309 , 1840.75183 , 1850.71033 , 1860.66858 , 1870.62659 , 1871.7843 , 1865.96375 , 1876.02527 , 1886.08459 , 1896.14148 , 1906.19604 , 1916.24841 , 1926.29846 , 1936.34607 , 1946.39148 , 1956.43457 , 1966.47534 , 1976.51379 , 1986.55005 , 1996.58386 , 2006.61536 , 2016.64465 , 2026.67163 , 2036.69617 , 2046.71863 , 2056.73877 , 2066.75635 , 2076.77173 , 2086.78491 , 2096.79565 , 2106.8042 , 2116.8103 , 2126.81421 , 2136.81567 , 2146.81494 , 2156.81201 , 2166.80664 , 2176.79907 , 2186.78906 , 2196.77661 , 2206.76221 , 2216.74536 , 2226.72607 , 2236.70459 , 2246.68066 , 2256.65454 , 2266.62622 , 2276.59546 , 2286.5625 , 2296.5271 , 2306.4895 , 2316.44946 , 2326.40723 , 2336.36255 , 2346.31567 , 2356.2666 , 2366.21509 , 2376.16113 , 2386.10522 , 2396.04663 , 2405.98608 , 2415.9231 , 2425.85767 , 2435.79004 , 2445.71997 , 2455.64771 , 2465.57324 , 2475.49634 , 2485.41724 , 2495.33569])

B) Using the indexes you found in part A, get those two bands from the reflectance file and assign them each to a variable (you can either read them from the file again with rasterio or you can get them using an index from the reflectance array). The ~857 band can be called nir and the ~645 band can be called red.

C) Use the formula from the question introduction to calculate NDVI

D) Save your new NDVI band to a file. Be sure to update the metadata to specify that there is now only 1 band in the file.