Merge pull request #64 from mxochicale/rasterio

adds examples of geospatial tools under `playground/geospatial` #63

.gitignore

@@ -18,7 +18,8 @@ checkpoints/
 ## ignores train_results
 **/train_results/*
 #!**\train_results\test1\eval.py
-## ignores MNIST data paths
+
+## ignores data paths
 **/data/**
 
 ## ignores vgg16 models

playground/geospatial/00_download_data.bash

@@ -0,0 +1,11 @@
+cd $HOME/repositories/mxochicale/code/playground/rasterio
+mkdir -p data && cd data
+
+#Download RBG.byte.tif and tif from googleapis
+wget https://github.com/rasterio/rasterio/raw/main/tests/data/RGB.byte.tif
+wget https://storage.googleapis.com/gcp-public-data-landsat/LC08/01/042/034/LC08_L1TP_042034_20170616_20170629_01_T1/LC08_L1TP_042034_20170616_20170629_01_T1_B4.TIF
+
+#Download spatial-vector-lidar data subset (~172 MB)
+wget https://ndownloader.figshare.com/files/12459464 -O spatial-vector-lidar.zip
+
+

playground/geospatial/01_raster_data.py

@@ -0,0 +1,75 @@
+#https://rasterio.readthedocs.io/en/stable/topics/plotting.html
+#https://github.com/rasterio/rasterio?tab=readme-ov-file#example
+
+import numpy as np
+import rasterio
+import matplotlib.pyplot as plt 
+from rasterio.plot import show
+from rasterio.plot import show_hist
+
+#Rasterio gives access to properties of a geospatial raster file.
+with rasterio.open('data/RGB.byte.tif') as src:
+    print(src.width, src.height)
+    print(src.crs)
+    print(src.transform)
+    print(src.count)
+    print(src.indexes)
+
+## Plotting three color bands and its histogram
+with rasterio.open('data/RGB.byte.tif') as src:
+    fig, (axrgb, axhist) = plt.subplots(1, 2, figsize=(14,7))
+    show(src, ax=axrgb)
+    show_hist(src, bins=50, histtype='stepfilled',
+          lw=0.0, stacked=False, alpha=0.3, ax=axhist)
+    plt.show()
+
+
+
+## Plotting three color bands
+with rasterio.open('data/RGB.byte.tif') as src:
+    fig, (axr, axg, axb) = plt.subplots(1,3, figsize=(21,7))
+    show((src, 1), ax=axr, cmap='Reds', title='red channel')
+    show((src, 2), ax=axg, cmap='Greens', title='green channel')
+    show((src, 3), ax=axb, cmap='Blues', title='blue channel')
+    plt.show()
+
+## Plotting countours
+with rasterio.open('data/RGB.byte.tif') as src:
+    fig, ax = plt.subplots(1, figsize=(12, 12))
+    show((src, 1), cmap='Greys_r', interpolation='none', ax=ax)
+    show((src, 1), contour=True, ax=ax)
+    plt.show()
+
+
+#A rasterio dataset also provides methods for getting read/write windows (like extended array slices) given georeferenced coordinates.
+with rasterio.open('data/RGB.byte.tif') as src:
+    window = src.window(*src.bounds)
+    print(window)
+    print(src.read(window=window).shape)
+
+# Read raster bands directly to Numpy arrays.
+with rasterio.open('data/RGB.byte.tif') as src:
+    r, g, b = src.read()
+    # Write the product as a raster band to a new 8-bit file. For
+    # the new file's profile, we start with the meta attributes of
+    # the source file, but then change the band count to 1, set the
+    # dtype to uint8, and specify LZW compression.
+    profile = src.profile
+    profile.update(dtype=rasterio.uint8, count=1, compress='lzw')
+
+# Combine arrays in place. Expecting that the sum will
+# temporarily exceed the 8-bit integer range, initialize it as
+# a 64-bit float (the numpy default) array. Adding other
+# arrays to it in-place converts those arrays "up" and
+# preserves the type of the total array.
+total = np.zeros(r.shape)
+#print(total.shape) #(718, 791)
+
+for band in r, g, b:
+    total += band
+
+total /= 3
+
+with rasterio.open('data/example-total.tif', 'w', **profile) as dst:
+    dst.write(total.astype(rasterio.uint8), 1)
+

playground/geospatial/02_raster_data.py

@@ -0,0 +1,46 @@
+#https://www.earthdatascience.org/workshops/gis-open-source-python/open-lidar-raster-python/
+
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import rasterio as rio
+from rasterio.plot import show
+from rasterio.plot import show_hist
+from rasterio.mask import mask
+import earthpy.plot as ep
+
+# define path to digital terrain model
+#sjer_dtm_path = "data/spatial-vector-lidar/california/neon-soap-site/2013/lidar/SOAP_lidarDTM.tif"
+sjer_dtm_path = "data/spatial-vector-lidar/california/neon-sjer-site/2013/lidar/SJER_lidarDSM.tif"
+
+# open raster data
+lidar_dem = rio.open(sjer_dtm_path)
+
+
+# read in all of the data without specifying a band
+with rio.open(sjer_dtm_path) as src:
+    print(src.bounds) # optional - view spatial extent
+    lidar_dem_im = src.read(masked= True) # convert / read the data into a numpy array:
+    #print(lidar_dem_im.shape) #(1, 1516, 3292)
+
+# specify a band so you get a 2 dimensional image array
+with rio.open(sjer_dtm_path) as src:
+    lidar_dem_im = src.read(1, masked= True) # convert / read the data into a numpy array:
+    sjer_ext = rio.plot.plotting_extent(src)
+    print(lidar_dem_im.shape) #(1516, 3292)
+    print(sjer_ext)
+
+ep.plot_bands(lidar_dem_im,
+              cmap='viridis_r',
+              extent=sjer_ext,
+              title="Lidar Digital Elevation Model \n Pre 2013 Boulder Flood | Lee Hill Road",
+              scale=False)
+plt.show()
+
+# create histogram of data
+ep.hist(lidar_dem_im,
+        bins=100,
+        title="Lee Hill Road - Digital Elevation (terrain) Model - \nDistribution of Elevation Values")
+plt.show()
+
+

playground/geospatial/03_shapefiles.py

@@ -0,0 +1,48 @@
+#https://www.earthdatascience.org/workshops/gis-open-source-python/intro-vector-data-python/
+
+import os
+import matplotlib.pyplot as plt
+import geopandas as gpd
+import earthpy as et
+
+sjer_plot_locations = gpd.read_file('data/spatial-vector-lidar/california/neon-sjer-site/vector_data/SJER_plot_centroids.shp')
+
+#Shapefile Structure
+#There are 3 key files associated with any and all shapefiles:
+#.shp: the file that contains the geometry for all features.
+#.shx: the file that indexes the geometry.
+#.dbf: the file that stores feature attributes in a tabular format.
+
+
+PRINT=True
+#PRINT=False
+if PRINT:
+    # view  the top 6 lines of attribute table of data
+    #print(sjer_plot_locations.head(6))
+    print(sjer_plot_locations)
+    #print(type(sjer_plot_locations)) #<class 'geopandas.geodataframe.GeoDataFrame'>
+
+    #view the spatial extent
+    #print(sjer_plot_locations.total_bounds) #[ 254738.618 4107527.074  258497.102 4112167.778]
+    #The spatial extent of a shapefile or geopandas GeoDataFrame represents the geographic "edge" or 
+    #location that is the furthest north, south east and west. 
+    #Thus is represents the overall geographic coverage of the spatial object. 
+    #Image Source: National Ecological Observatory Network (NEON)    
+    #print(sjer_plot_locations.crs) #EPSG:32611
+    #print(sjer_plot_locations.geom_type)
+    #print(sjer_plot_locations.shape) #(18, 6)
+
+
+## quickly plot the data adding a legend
+fig, ax = plt.subplots(figsize = (10,10))
+sjer_plot_locations.plot(column='plot_type', 
+                         categorical=True, 
+                         legend=True, 
+                         figsize=(10,6),
+                         marker='*', 
+                         markersize=45, 
+                         cmap="Set2", ax=ax);
+ax.set_title('SJER Plot Locations\nMadera County, CA', fontsize=16);
+plt.show()
+
+

playground/geospatial/04_cropping_raster_data.py

@@ -0,0 +1,80 @@
+#https://www.earthdatascience.org/workshops/gis-open-source-python/crop-raster-data-in-python/
+
+import os
+import numpy as np
+import rasterio as rio
+from rasterio.plot import show
+from rasterio.mask import mask
+import matplotlib.pyplot as plt
+import earthpy as et
+import earthpy.plot as ep
+import earthpy.spatial as es
+import geopandas as gpd
+
+soap_chm_path = 'data/spatial-vector-lidar/california/neon-soap-site/2013/lidar/SOAP_lidarCHM.tif'
+
+# open the lidar chm
+with rio.open(soap_chm_path) as src:
+    lidar_chm_im = src.read(masked=True)[0]
+    extent = rio.plot.plotting_extent(src)
+    soap_profile = src.profile
+
+#print(type(lidar_chm_im))
+#<class 'numpy.ma.core.MaskedArray'>
+
+#ep.plot_bands(lidar_chm_im,
+#               cmap='terrain',
+#               extent=extent,
+#               title="Lidar Canopy Height Model (CHM)\n NEON SOAP Field Site",
+#               cbar=False);
+#
+
+#Open Vector Layer
+crop_extent_soap = gpd.read_file('data/spatial-vector-lidar/california/neon-soap-site/vector_data/SOAP_crop2.shp')
+
+#explore the coordinate reference system (CRS)
+print('crop extent crs: ', crop_extent_soap.crs)
+print('lidar crs: ', soap_profile['crs'])
+
+# plot the data
+fig, ax = plt.subplots(figsize = (6, 6))
+crop_extent_soap.plot(ax=ax)
+ax.set_title("Shapefile Imported into Python \nCrop Extent for Soaproot Saddle Field Site", 
+             fontsize = 16)
+ax.set_axis_off()
+plt.show()
+
+fig, ax = plt.subplots(figsize=(10, 10))
+ep.plot_bands(lidar_chm_im,
+              cmap='terrain',
+              extent=extent,
+              ax=ax,
+              cbar=False)
+crop_extent_soap.plot(ax=ax, alpha=.6, color='g');
+plt.show()
+
+#To crop the data,use the crop_image function in earthpy.spatial.
+with rio.open(soap_chm_path) as src:
+    lidar_chm_crop, soap_lidar_meta = es.crop_image(src, crop_extent_soap)
+
+# Update the metadata to have the new shape (x and y and affine information)
+soap_lidar_meta.update({"driver": "GTiff",
+                 "height": lidar_chm_crop.shape[0],
+                 "width": lidar_chm_crop.shape[1],
+                 "transform": soap_lidar_meta["transform"]})
+
+# generate an extent for the newly cropped object for plotting
+cr_ext = rio.transform.array_bounds(soap_lidar_meta['height'], 
+                                            soap_lidar_meta['width'], 
+                                            soap_lidar_meta['transform'])
+
+bound_order = [0,2,1,3]
+cr_extent = [cr_ext[b] for b in bound_order]
+cr_extent, crop_extent_soap.total_bounds
+
+# mask the nodata and plot the newly cropped raster layer
+lidar_chm_crop_ma = np.ma.masked_equal(lidar_chm_crop[0], -9999.0) 
+# plot
+ep.plot_bands(lidar_chm_crop_ma, cmap='terrain', cbar=False)
+
+

playground/geospatial/README.md

@@ -0,0 +1,21 @@
+# Geospatial Analysis Tools (rasterio, earthpy, etc)
+
+## Download data
+```
+cd $HOME/repositories/mxochicale/code/playground/geospatial
+bash 00_download_data.bash
+```
+
+## Scripts
+```
+mamba activate eVE
+python *.py
+```
+
+## References
+https://geohackweek.github.io/raster/04-workingwithrasters/
+https://www.tomasbeuzen.com/python-for-geospatial-analysis/chapters/chapter1_intro-to-spatial.html
+https://apoorvalal.github.io/notebook/r-py-notes/cu_earthlab/cu_earthlab_intermediate.html 
+https://earthpy.readthedocs.io/en/latest/index.html
+
+

pyVEs/eVE.yml

@@ -27,6 +27,9 @@ dependencies:
      - tqdm #https://github.com/tqdm/tqdm/tags
      - pylint #https://github.com/pylint-dev/pylint/tags
      - seaborn #https://github.com/mwaskom/seaborn/tags
+     - rasterio
+     - earthpy #https://github.com/earthlab/earthpy
+     - geopandas #https://github.com/geopandas/geopandas
      ### VERSIONS of pyside6: https://pypi.org/project/PySide6/#history
      #- pyside6>=6.4.2
      ### VERSIONS of VTK https://gitlab.kitware.com/vtk/vtk/-/tags