Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
###############################################################################
# Import Example Data
# -------------------
#
# To get started, make sure your directory is set. Then, create a stack from all
# of the Landsat .tif files (one per band).
# Get data for example
data = et.data.get_data("vignette-landsat")
# Set working directory
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
# Stack the Landsat 8 bands
# This creates a numpy array with each "layer" representing a single band
# You can use the nodata parameter to mask nodata values
landsat_path = glob(os.path.join("data", "vignette-landsat", "*band*.tif"))
landsat_path.sort()
array_stack, meta_data = es.stack(landsat_path, nodata=-9999)
###############################################################################
# Plot All Histograms in a Stack With Custom Titles and Colors
# -------------------------------------------------------------
#
# You can create histograms for each band with unique colors and titles by
# first creating a list of colors and titles that will be provided to the
# ep.hist() function. The list for titles must contain the same number of
# strings as there are bands in the stack. You can also modify the colors for
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
import earthpy.mask as em
# Get data and set your home working directory
data = et.data.get_data("cold-springs-fire")
###############################################################################
# Import Example Data
# ------------------------------
# To get started, make sure your directory is set. Create a stack from all of the
# Landsat .tif files (one per band) and import the ``landsat_qa`` layer which provides
# the locations of cloudy and shadowed pixels in the scene.
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
# Stack the landsat bands
# This creates a numpy array with each "layer" representing a single band
landsat_paths_pre = glob(
"data/cold-springs-fire/landsat_collect/LC080340322016070701T1-SC20180214145604/crop/*band*.tif"
)
landsat_paths_pre.sort()
arr_st, meta = es.stack(landsat_paths_pre)
# Import the landsat qa layer
with rio.open(
"data/cold-springs-fire/landsat_collect/LC080340322016070701T1-SC20180214145604/crop/LC08_L1TP_034032_20160707_20170221_01_T1_pixel_qa_crop.tif"
) as landsat_pre_cl:
landsat_qa = landsat_pre_cl.read(1)
landsat_ext = plotting_extent(landsat_pre_cl)
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
import rasterio as rio
# Download the data needed for this vignette
data = et.data.get_data("vignette-elevation")
####################################################################################
# Open up the DEM
# ---------------
# To begin, open your DEM layer as a numpy array using Rasterio. Below you set all
# terrain values < 0 to ``nan``. Next, plot the data using ``ep.plot_bands()``.
# Set the home directory and get the data for the exercise
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
dtm = "data/vignette-elevation/pre_DTM.tif"
# Open the DEM with Rasterio
with rio.open(dtm) as src:
elevation = src.read(1)
# Set masked values to np.nan
elevation[elevation < 0] = np.nan
# Plot the data
ep.plot_bands(
elevation,
scale=False,
cmap="gist_earth",
title="DTM Without Hillshade",
figsize=(10, 6),
)
data = et.data.get_data("spatial-vector-lidar")
###############################################################################
# Open Files with GeoPandas and Reproject the Data
# -------------------------------------------------
#
# Start by setting your working directory. Then, import the data files to
# GeoDataFrames using GeoPandas.
#
# Recall that the data must be in the same CRS in order to use the
# ``clip_shp()`` function. If the data are not in the same CRS, be sure to use
# the ``to_crs()`` function from GeoPandas to match the projects between the
# two objects, as shown below.
# Set your home environment
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
# Open both files with GeoPandas
road_path = os.path.join(
"data",
"spatial-vector-lidar",
"global",
"ne_10m_roads",
"ne_10m_n_america_roads.shp",
)
roads = gpd.read_file(road_path)
country_path = os.path.join(
"data", "spatial-vector-lidar", "usa", "usa-boundary-dissolved.shp"
)
country_boundary = gpd.read_file(country_path)
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
# Get data and set your home working directory
data = et.data.get_data("cold-springs-fire")
###############################################################################
# Import Example Data
# -------------------
#
# To get started, make sure your directory is set. Create a stack from all of the
# Landsat .tif files (one per band).
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
# Stack the Landsat 8 bands
# This creates a numpy array with each "layer" representing a single band
landsat_path = glob(
"data/cold-springs-fire/landsat_collect/LC080340322016072301T1-SC20180214145802/crop/*band*.tif"
)
landsat_path.sort()
arr_st, meta = es.stack(landsat_path)
###############################################################################
# Calculate Normalized Difference Vegetation Index (NDVI)
# -------------------------------------------------------
# You can calculate NDVI for your dataset using the
# ``normalized_diff`` function from the ``earthpy.spatial`` module.
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
# Get data and set your home working directory
data = et.data.get_data("cold-springs-fire")
###############################################################################
# Import Example Data
# -------------------
#
# To get started, make sure your directory is set. Then, create a stack from all of
# the Landsat .tif files (one per band).
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
# Stack the Landsat 8 bands
# This creates a numpy array with each "layer" representing a single band
landsat_path = glob(
"data/cold-springs-fire/landsat_collect/LC080340322016072301T1-SC20180214145802/crop/*band*.tif"
)
landsat_path.sort()
arr_st, meta = es.stack(landsat_path)
###############################################################################
# Calculate Normalized Difference Vegetation Index (NDVI)
# -------------------------------------------------------
#
# You can calculate NDVI for your dataset using the
# ``normalized_diff`` function from the ``earthpy.spatial`` module.
# into one file. With EarthPy, you can create an image stack from all of the
# Landsat .tif files (one per band) in a directory using the ``stack()`` function
# from the ``earthpy.spatial`` module.
###################################################################################
# Error found on Windows systems
# -------------------------------
# .. note::
# If you are running this script on a Windows system, there is a
# known bug with ``.to_crs()``, which is used in this script. If an error
# occurs, you have to reset your os environment with the command
# ``os.environ["PROJ_LIB"] = r"path-to-share-folder-in-environment"``.
# Get sample data from EarthPy and set working directory
data_path = et.data.get_data("vignette-landsat")
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
# Get list of bands and sort by ascending band number
landsat_bands_data_path = "data/vignette-landsat/LC08_L1TP_034032_20160621_20170221_01_T1_sr_band*[1-7]*.tif"
stack_band_paths = glob(landsat_bands_data_path)
stack_band_paths.sort()
# Create image stack and apply nodata value for Landsat
arr_st, meta = es.stack(stack_band_paths, nodata=-9999)
###############################################################################
# Plot RGB Composite Image
# --------------------------
# You can use the ``plot_rgb()`` function from the ``earthpy.plot`` module to quickly
# plot three band composite images. For RGB composite images, you will plot the red,
# green, and blue bands, which are bands 4, 3, and 2, respectively, in the image
# stack you created. Python uses a zero-based index system, so you need to subtract
plt.show()
#############################################################################
# Plot Boundary Over Composite Image
# -----------------------------------
# .. note::
# If you are on windows, you may need to add the crs issue discussed above
# here!
#
# You can overlay a polygon boundary on top of an RGB plot created with EarthPy.
# To begin, the raster data and the boundary need to be in the same
# Coordinate Reference System (CRS). You can reproject the boundary layer to
# match the CRS of the image by getting the CRS of the image from the Rasterio
# profile object and passing that CRS to the ``to_crs`` method from GeoPandas.
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
# Open polygon boundary using GeoPandas
bound = gpd.read_file(
"data/vignette-landsat/vector_layers/fire-boundary-geomac/co_cold_springs_20160711_2200_dd83.shp"
)
# Reproject boundary to match CRS of the Landsat images
with rio.open(stack_band_paths[0]) as raster_crs:
raster_profile = raster_crs.profile
bound_utm13N = bound.to_crs(raster_profile["crs"])
################################################################################
# Create a Plot With the Boundary overlayed on the RGB Image
# ----------------------------------------------------------
# You can plot a polygon boundary over an image by creating a raster extent
# for the plot using the ``plotting_extent`` function from ``rasterio.plot``.
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
###############################################################################
# Import Example Data
# -------------------
#
# To get started, make sure your directory is set. Then, create a stack from all of
# the Landsat .tif files (one per band).
# Get data for example
data = et.data.get_data("vignette-landsat")
# Set working directory
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
# Stack the Landsat 8 bands
# This creates a numpy array with each "layer" representing a single band
# You can use the nodata= parameter to mask nodata values
landsat_path = glob(
"data/vignette-landsat/LC08_L1TP_034032_20160621_20170221_01_T1_sr_band*_crop.tif"
)
landsat_path.sort()
array_stack, meta_data = es.stack(landsat_path, nodata=-9999)
###############################################################################
# Plot All Bands in a Stack
# --------------------------
#
# When you give ``ep.plot_bands()`` a three dimensional numpy array,
# it will plot all layers in the numpy array. You can create unique titles for
os.mkdir(output_dir)
raster_out_path = os.path.join(output_dir, "raster.tiff")
####################################################################################
# Stack the Bands
# ---------------------------
# The stack function has an optional output argument, where you can write the raster
# to a tiff file in a folder. If you want to use this functionality, make sure there
# is a folder to write your tiff file to.
# The Stack function also returns two object, an array and a RasterIO profile. Make
# sure to be catch both in variables.
# Stack Landsat bands
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
array, raster_prof = es.stack(stack_band_paths, out_path=raster_out_path)
####################################################################################
# Create Extent Object
# --------------------------------
# To get the raster extent, use the ``plotting_extent`` function on the
# array from ``es.stack()`` and the Rasterio profile or metadata object. The function
# needs a single
# layer of a numpy array, which is why we use ``arr[0]``. The function also
# needs the spatial transformation for the Rasterio object, which can be acquired by accessing
# the ``"transform"`` key within the Rasterio Profile.
extent = plotting_extent(array[0], raster_prof["transform"])
################################################################################
# Plot Un-cropped Data