diff --git a/docs/exercises/raster_and_vector.html b/docs/exercises/raster_and_vector.html index 728f5cf..959b27b 100644 --- a/docs/exercises/raster_and_vector.html +++ b/docs/exercises/raster_and_vector.html @@ -881,7 +881,7 @@

86  +
# Import modules
 import numpy as np
 import geopandas as gpd
@@ -891,12 +891,12 @@ 
86  
 
Read vector layer

-

+

 
# Read file with watersheds
 vector_path = '../datasets/spatial/konza_watersheds.geojson'
 gdf = gpd.read_file(vector_path)
 gdf.head(3)

-

+

 

 
 
@@ -958,7 +958,7 @@ 
Read vector layer

 

 

 

-

+

 
# Display coordiante reference system of vector layer
 print(gdf.crs)

 

@@ -967,11 +967,11 @@ 
Read vector layer

 

 
The first step to perform spatial operations between vector and raster layers is to have both data sources in the same coordinate reference system. Our vector layer has geometries in projected coordinates (meters), but the raster datasets will be in geographic coordinates.

 
To perform operations that involve distance computations is best to have both datasets in projected coordiantes. To perform other operations like clipping layers, then both geographic and projected coordinates should perform similarly. In this example we will convert the vector layer to geographic coordinates.

-

+

 
# Change coordinate reference system from projected to geographic
 gdf.to_crs(epsg=4326, inplace=True)

 

-

+

 
# Visualize vector map of watersheds
 gdf.plot(facecolor='None', edgecolor='k')
 plt.title('Konza Watersheds')
@@ -992,11 +992,11 @@ 
Read raster layer

 
  • Band 3: Sentinel B08 (near infrared, 842 nm)
    • 
 
 
    The image was downloaded using Sentinel Hub: https://www.sentinel-hub.com/explore/eobrowser/
    
-
    
+
    
 
    # Load raster layer
 raster = xr.open_dataarray('../datasets/spatial/2024-04-03_Sentinel-2_L2A_B03_B04_B08-16bit.tiff')
 raster
    
-
    
+
    
 
 
    
 
@@ -1372,8 +1372,8 @@ 
    Read raster layer
    
     AREA_OR_POINT:           Area
     TIFFTAG_RESOLUTIONUNIT:  1 (unitless)
     TIFFTAG_XRESOLUTION:     1
-    TIFFTAG_YRESOLUTION:     1
  • AREA_OR_POINT :
    Area
    TIFFTAG_RESOLUTIONUNIT :
    1 (unitless)
    TIFFTAG_XRESOLUTION :
    1
    TIFFTAG_YRESOLUTION :
    1
    • -
    +
    # Display CRS for raster layer
 print(raster.rio.crs)
    EPSG:4326
    -
    +
    # Create separate variables to easily keep track of each band
-B4 = raster[1] # Red band
-B8 = raster[2] # NIR band
    +red = raster[1] # Red band +nir = raster[2] # NIR band
    -
    +
    # Visualize raster files (1 band each)
 plt.figure(figsize=(12,4))
 
 plt.subplot(1,2,1)
-plt.title('Band 4')
-B4.plot(cmap='Spectral')
+plt.title('Red band')
+red.plot(cmap='Spectral')
 
 plt.subplot(1,2,2)
-plt.title('Band 8')
-B8.plot(cmap='Spectral')
+plt.title('NIR band')
+nir.plot(cmap='Spectral')
 
 plt.show()
    @@ -1427,28 +1427,28 @@

    Read raster layer

    Compute NDVI

    Normalized Difference Vegetation Index (NDVI) is a measure typically used in remote sensing to quantify and assess vegetation health and density. It is calculated using near-infrared (NIR) and red light reflectance values over an area as follows:

    -

    NDVI = \frac{Red - NIR}{Red+NIR}

    +

    NDVI = \frac{NIR-Red}{NIR+Red}

    NDVI values range from -1 to 1, where higher values indicate healthier and denser vegetation, typically falling between 0.2 and 0.8 for vegetation. Bare soil and rocks fall have values ranging from 0 to 0.2. Bodies of water like ponds, likes, rivers, and oceans typically have negative values close to -1.

    -
    +
    # Calculate NDVI
-NDVI = (B8 - B4)/(B8 + B4)
    +ndvi = (nir - red)/(nir + red)

    Now that we have the NDVI raster layer, here a few additional commands that could be useful:

    # Get NDVI values as Numpy array
-NDVI.values
+ndvi.values
 
 # Get latitude values as Numpy array
-NDVI.coords['y'].values
+ndvi.coords['y'].values
 
 # Get a 3 by 3 slice
-NDVI[0:3,0:3]
+ndvi[0:3,0:3]
 
 # Find NDVI value for nearest point
-NDVI.sel(x=-96.664356, y=39.15406606, method='nearest')
    +ndvi.sel(x=-96.664356, y=39.15406606, method='nearest')

    Define custom colormap

    -
    +
    # Palette of colors NDVI
 hex_palette = ['#CE7E45', '#DF923D', '#F1B555', '#FCD163', '#99B718', '#74A901',
              '#66A000', '#529400', '#3E8601', '#207401', '#056201', '#004C00', '#023B01',
@@ -1459,17 +1459,17 @@ 
    Define custom color
 ndvi_cmap.set_bad('#FEFEFE')
 ndvi_cmap.set_under('#0000FF')
 ndvi_cmap
    -
    +
    from_list
    from_list colormap
    under
    bad
    over
    -
    +
    # Add new colormap to list of Matplotlib's colormaps
 plt.colormaps.register(cmap=ndvi_cmap, name='ndvi')
    -
    +
    # Visualize NDVI for the entire area of the image
-NDVI.plot(cmap='ndvi', vmin=0.0, vmax=1.0, 
+ndvi.plot(cmap='ndvi', vmin=0.0, vmax=1.0, 
           cbar_kwargs={'label':'NDVI', 'shrink':0.5})
 plt.title('NDVI')
 plt.xlabel('Longitude')
@@ -1482,18 +1482,18 @@ 
    Define custom color

    Clip NDVI raster to watersheds

    -
    +
    # Define a function to clip the raster with each polygon and return a numpy array
 clip_fn = lambda polygon, R: R.rio.clip([polygon.geometry], 
                                         crs=R.rio.crs, 
                                         all_touched=True)
 
 # Apply the function to each row in the GeoDataFrame to create a new 'clipped_raster' column
-gdf['clipped_raster'] = gdf.apply(lambda row: clip_fn(row, NDVI), axis=1)
+gdf['clipped_raster'] = gdf.apply(lambda row: clip_fn(row, ndvi), axis=1)
 
 # Inspect resulting GeoDataframe
 gdf.head(3)
    -
    +
    @@ -1563,7 +1563,7 @@

    Clip NDVI r

    Visualize a specific clipped watershed

    We will choose the K1B watershed to show the NDVI for the entire area.

    -
    +
    # Select watershed
 idx = gdf['NAME_1'] == 'K1B'
 row = gdf[idx].index[0]
@@ -1572,7 +1572,7 @@ 
    Vis
 
    31
    -
    +
    # Create figure of selected watershed
 fig, ax = plt.subplots(figsize=(6, 6))
 gdf.loc[ [row], 'geometry'].boundary.plot(ax=ax, edgecolor='k')
@@ -1592,7 +1592,7 @@ 
    Vis

    Compute mean NDVI for each watershed

    -
    +
    # Create empty list to append watershed NDVI values
 ndvi_avg = []
 
@@ -1605,7 +1605,7 @@ 
    Compu
 
 # Inspect results
 gdf.head(3)
    -
    +
    @@ -1678,7 +1678,7 @@

    Compu

    Create static map

    -
    +
    # Create figure with mean NDVI values for each watershed
 gdf.plot(column='ndvi_avg', edgecolor='k', cmap='ndvi', 
          legend=True, vmin=0.0, vmax=1.0, legend_kwds={'label':'NDVI', 'shrink':0.5})
@@ -1694,12 +1694,12 @@ 
    Create static map

    Create interactive map

    -
    +
    # Create interactive map
 gdf.iloc[:,:-1].explore(column='ndvi_avg', k=10, cmap='ndvi',
                         style_kwds={'stroke':True,'color':'black','width':1},
                         legend_kwds={'caption':'NDVI'})
    -
    +
    Make this Notebook Trusted to load map: File -> Trust Notebook
    " ], "text/plain": [ - "" + "" ] }, - "execution_count": 20, + "execution_count": 43, "metadata": {}, "output_type": "execute_result" }