From ab486b92a0ffbdd78ca58088ef963cd1b3f731ff Mon Sep 17 00:00:00 2001 From: mxochicale Date: Sat, 16 Mar 2024 15:46:32 +0000 Subject: [PATCH] adds examples of geospatial tools under `playground/geospatial` #63 --- .gitignore | 3 +- playground/geospatial/00_download_data.bash | 11 +++ playground/geospatial/01_raster_data.py | 75 +++++++++++++++++ playground/geospatial/02_raster_data.py | 46 +++++++++++ playground/geospatial/03_shapefiles.py | 48 +++++++++++ .../geospatial/04_cropping_raster_data.py | 80 +++++++++++++++++++ playground/geospatial/README.md | 21 +++++ pyVEs/eVE.yml | 3 + 8 files changed, 286 insertions(+), 1 deletion(-) create mode 100644 playground/geospatial/00_download_data.bash create mode 100644 playground/geospatial/01_raster_data.py create mode 100644 playground/geospatial/02_raster_data.py create mode 100644 playground/geospatial/03_shapefiles.py create mode 100644 playground/geospatial/04_cropping_raster_data.py create mode 100644 playground/geospatial/README.md diff --git a/.gitignore b/.gitignore index 6f282fa..c7bb216 100644 --- a/.gitignore +++ b/.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 diff --git a/playground/geospatial/00_download_data.bash b/playground/geospatial/00_download_data.bash new file mode 100644 index 0000000..e668aa9 --- /dev/null +++ b/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 + + diff --git a/playground/geospatial/01_raster_data.py b/playground/geospatial/01_raster_data.py new file mode 100644 index 0000000..6e58e10 --- /dev/null +++ b/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) + diff --git a/playground/geospatial/02_raster_data.py b/playground/geospatial/02_raster_data.py new file mode 100644 index 0000000..0cbae36 --- /dev/null +++ b/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() + + diff --git a/playground/geospatial/03_shapefiles.py b/playground/geospatial/03_shapefiles.py new file mode 100644 index 0000000..11656a1 --- /dev/null +++ b/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)) # + + #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() + + diff --git a/playground/geospatial/04_cropping_raster_data.py b/playground/geospatial/04_cropping_raster_data.py new file mode 100644 index 0000000..55509c8 --- /dev/null +++ b/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)) +# + +#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) + + diff --git a/playground/geospatial/README.md b/playground/geospatial/README.md new file mode 100644 index 0000000..9894667 --- /dev/null +++ b/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 + + diff --git a/pyVEs/eVE.yml b/pyVEs/eVE.yml index 11a1398..15897a0 100644 --- a/pyVEs/eVE.yml +++ b/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