After completing this tutorial, you will be able to:
- Find and download data from the USGS Earth Explorer Website
- Filter data by cloud cover to find datasets with the least amount of clouds for a study area.
What you need
You will need a computer with internet access to complete this lesson and the data for week 6 / 7 of the course.
In class this week, we will review how to grab data from the Earth Explorer website. The Earth Explorer website is a data portal run by the USGS. Here you can find many different types of remote sensing and other data for both the US and in some cases, the globe.
IMPORTANT: Be sure to order your data several days ahead of time or else you won’t have it in time to finish this assignment.
The Steps: Earth Explorer Data Download
Define study area (AOI)
When searching for data, the first thing we need to do is to define our area of interest (AOI). Our AOI is defined by the boundary of the fire extent. We could type in the x,y vertices of each corner of the boundary, but if we have an Earth Explorer account, we can upload a ZIPPED up shapefile that contains the boundary instead!
Important: Be sure to use a square extent. If you have too many vertices in your extent polygon, the website won’t accept it as an extent file.
- Zip up extent file that you want to use. Be sure to use a square extent, if you have too many vertices it won’t work. Lucky for us there is a zip file already zipped up and ready to go in our week6 data!
- Next, go to the Earth Explorer website. Login. If you don’t have a login already, create an account.
Now, it’s time to search for data.
- In the search criteria, click on shapefile tab. Select the zip file above as the shapefile that represents the SPATIAL EXTENT of our study area - the cold springs fire site.
- At the bottom of the search criteria window, select a range of dates. A month before and after the fire is a nice starting point.
- Next click on the Data sets tab. Notice that there are a lot of different data available from Earth Explorer! We are interested in Landsat - specifically Landsat 8. You can find Landsat in the Landsat archive drop down. Expand that drop down to find:
- Pre-Collection -> Landsat Surface Reflectance - L8 OLI/TIRS
- Next select the Additional Criteria tab. Here is where you can limit results by % cloud cover. Let’s start with Less than 20% cloud cover and see what we get as data results.
- Finally click on the Results tab. Here you see all of the scenes available for “order” from the website that cover our study area.
- Notice that you can click on the icons below the scene to see the scene itself rendered on the map and to see the footprint (or extent) of the scene relative to our study area.
- Pick a scene that is
- closest to the pre-fire date (July 10 2016) and also that has the least amount of cloud cover close to our study area.
Order your data
- Click the shopping cart icon to add the data to your cart.
- Click on “item basket” in the upper right hand corner of your browser to see what you have ordered.
- Click on Proceed to Checkout
- Then finally, click on Submit Order
IMPORTANT: It will take a few days for the link that you can use to download your data to be emailed to your account. Order now!
In this case, I downloaded a scene very close to Julian day 189.
Import new scene
First, let’s import our new data and create a raster stack. The code is hidden because you already know how to do this!
## Error in x[]: subscript out of bounds ## Error in plotRGB(all_landsat_bands_173_st, r = 4, g = 3, b = 2, stretch = "lin", : object 'all_landsat_bands_173_st' not found ## Error in box(col = "white"): plot.new has not been called yet
Next I plotted the fire boundary extent on top of my landsat image.
## Error in ogrListLayers(dsn = dsn): Cannot open data source ## Error in crs(all_landsat_bands_st): object 'all_landsat_bands_st' not found
## Error in plotRGB(all_landsat_bands_173_st, r = 4, g = 3, b = 2, stretch = "lin", : object 'all_landsat_bands_173_st' not found ## Error in box(col = "white"): plot.new has not been called yet ## Error in polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col, : plot.new has not been called yet
It’s hard to see but can you see the tiny YELLOW outline of our study area? This landsat scene is MUCH larger than our study area. We have 2 options
- Crop the data: this will make it easier to work with as it will be smaller. A good move.
- Plot only the study area extent: this is ok if we just want to plot our data and don’t need to do any additional processing on it.
Below i’ve plotted the cloud mask for the data that I downloaded. It looks like the data in our study area are cloud free. How do I know that?
## Error in crop(all_landsat_bands_173_st, fire_boundary_utm): object 'all_landsat_bands_173_st' not found ## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. (file does not exist) ## Error in crop(cloud_mask_173, fire_boundary_utm): object 'cloud_mask_173' not found ## Error in plot(cloud_mask_173_crop, main = "Cropped cloud mask layer for new Landsat scene", : object 'cloud_mask_173_crop' not found ## Error in polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col, : plot.new has not been called yet ## Error in legend([email protected]@xmax, [email protected]@ymax, : object 'cloud_mask_173_crop' not found
barplot(cloud_mask_173_crop, main = "cloud mask values \n all 0's") ## Error in barplot(cloud_mask_173_crop, main = "cloud mask values \n all 0's"): object 'cloud_mask_173_crop' not found
Given our data are all 0’s we can assume we downloaded the right scene! There are no clouds in our study area image. This means we don’t have to worry about masking.
# turn axes to white par(col.axis="white", col.lab="white", tck=0) # plot RGB plotRGB(all_landsat_bands_173_st, r=4, g=3, b=2, stretch="lin", main = "Final landsat scene with the fire extent overlayed", axes=T) ## Error in plotRGB(all_landsat_bands_173_st, r = 4, g = 3, b = 2, stretch = "lin", : object 'all_landsat_bands_173_st' not found box(col = "white") ## Error in box(col = "white"): plot.new has not been called yet plot(fire_boundary_utm, add=T, border="yellow") ## Error in polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col, : plot.new has not been called yet
Now we can proceed to calculate NBR on the pre-fire landsat image. How does it look?