In this exercise, you will learn how to read and write following types of spatial data in Python.
We will use following data set and could be found here.
- New York State county shape file (Polygon: NY_County_GCS)
- Soil Carbon Data (SOC) from Colorado, Kansas, New Mexico, and Wyoming (Point: GP_SOC_GCS.shp)
- Road network maps of Onondaga county, New York State (polyline: "Onondaga_Street_GCS.shp)
- 90 m SRTM DEM of Onondaga county, New York State (raster: Onondaga_DEM.tif)
- Landsat8 multibands Images
Set working directory¶
import os
path= "C:/Users/zahmed2/Documents/spatial_data_processing_and_analysis_python/Lesson_02_read_and_write_spatila_data_in_python/"
os.chdir(path)
print("Current Working Directory " , os.getcwd())
Vector Data¶
Vector data consists of "geometry" or "shape" of the locations which describe information of spatial objects on earth. There are three types vector data such as polygons, points and line or polylines which structures of the geometry consists of sets of coordinate pairs (x, y).
The shapefile, a interchangeably data format which are regulated by ESRI. It is one of the most common form of geospatial vector data used in GIS software and analyses. Three unique files are required for a shapefile, including: .shx (shape index format; this tags the shapefile with a position, so users can move it forward and backward among layers, a .shp (shape format, which stores geometric data, and a .dbf (attribute format file; which holds attributes (information) for the shapes in the file).
Polygon data¶
Polygons are two-dimensional geographical features covering a portion of the earth's surface, for example, forests, lakes, administrative boundaries, farmers' fields, or any other organizational unit the user defines. Polygons are important because their area and perimeter can be measured. Following example will how read ESRI polygon shape with four python libraries. We will use New York Sate County shape file (NY_County_GCS.shp).
Using Geopandas¶
import geopandas as gpd
# Set filepath
fp_poly_gcs="Data_02/NY_County_GCS.shp"
# read polygoan files with Geopandas
poly_gcs= gpd.GeoDataFrame.from_file(fp_poly_gcs)
# Check file type
type(poly_gcs)
# See data frame
poly_gcs.head()
Using PyShp¶
import shapefile as shp
shape = shp.Reader(fp_poly_gcs)
Using Fiona¶
import fiona
shape = fiona.open("Data_02/NY_County_GCS.shp")
print (shape.schema)
#first feature of the shapefile
first = shape.next()
print (first)
Osgeo- GDAL¶
from osgeo import ogr
file = ogr.Open(fp_poly_gcs)
shape = file.GetLayer(0)
print(shape)
Visiualization of Vector Map¶
#%matplotlib inline
import os
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import box
from mpl_toolkits.axes_grid1 import make_axes_locatable
# Import NY Ctate County Shape file
poly_gcs.plot()
### Change color of map
from matplotlib.colors import ListedColormap
cmap = ListedColormap(['grey'])
poly_gcs.plot(cmap=cmap)
## Plot catgorical data
fig, ax = plt.subplots(figsize=(8, 8))
poly_gcs.plot(column='COUNTY',
categorical=True,
legend=True,
ax=ax)
# adjust legend location
leg = ax.get_legend()
leg.set_bbox_to_anchor((1.25,1.25))
ax.set_axis_off()
plt.show()
## Plot Continues
# set a variable that will call whatever column we want to visualise on the map
variable = 'AREA'
# create figure and axes for Matplotlib
fig, ax = plt.subplots(1, figsize=(10, 6))
# create map
poly_gcs.plot(column=variable,
cmap='Blues',
linewidth=0.8,
ax=ax,
edgecolor='0.8')
# remove the axis
ax.axis('off')
# add a title
ax.set_title('NYS County Area')
# create an annotation for the data source
ax.annotate('Source: Cornell GIS, 2014',
xy=(0.1, .08),
xycoords='figure fraction',
horizontalalignment='left',
verticalalignment='top',
fontsize=12,
color='#555555')
# Create colorbar as a legend
sm = plt.cm.ScalarMappable(cmap='Blues')
# empty array for the data range
sm._A = []
# add the colorbar to the figure
cbar = fig.colorbar(sm)
Lines or Polylines¶
One-dimensional lines, also called polylines, are used to represent geographical features like rivers, roads, railroads, trails, and topographic lines. Note that these features are linear in nature and do not have area like polygons. Hence, they can measure distance. In this exercise we use road network maps of Onondaga county, New York State.
# Load road shape file
fp_road_gcs="Data_02/Ononda_Street_GCS.shp"
road_gcs= gpd.read_file(fp_road_gcs)
road_gcs.plot()
# See data frame
road_gcs.head()
## Plot road
fig, ax = plt.subplots(figsize=(14, 6))
road_gcs.plot(column='Types',
categorical=True,
legend=True,
ax=ax)
# adjust legend location
leg = ax.get_legend()
leg.set_bbox_to_anchor((1.5,1.05))
ax.set_axis_off()
plt.show()
Point data¶
Zero-dimensional points are important for geographical features like wells, soil-sampling locations that can be best expressed by a single reference point. Often point vector data consist off several points into a multi-point structure, with a single attribute record. For example, all the soil sampling points could be considered as a single geometry. The main drawback of a point feature, however, are that they cannot be used to make measurements (as you can with a polygon).
We used 650 soil samples from Colorado, Kansas, New Mexico, and Wyoming. These samples were collected by the United States Geological Survey (USGS) as a part of the USGS Geo-chemical Landscapes Project [Smith et al., 2011]. We will use shapefile() or readOGR() to read this vector data.
# Load point shape file
fp_point_gcs="Data_02/GP_Point_GCS.shp"
point_gcs= gpd.read_file(fp_point_gcs)
point_gcs.plot()
point_gcs.head()
Write Vector Data¶
GeoDataFrames can be exported to many different standard formats using the GeoDataFrame.to_file() method. For a full list of supported formats, type import fiona; fiona.supported_drivers.
Writing to Shapefile:¶
# poly_gcs.to_file(Data_02/NY_County.shp")
Writing to GeoJSON:¶
# poly_gcs.to_file("Data_02/NY_County_geo", driver='GeoJSON')
Raster Data¶
Unlike vector data, a raster data consists of cells or pixels organized into rows and columns as a matrix where each cell contains a value representing geographical feature on the earth. The size of the area on the surface that each pixel covers is known as the spatial resolution of the image. For instance, an image that has a 1 m spatial resolution means that each pixel in the image represents a 1 m x 1 m area. There are two types of raster data: continuous and discrete. An example of discrete raster data is Land use raster. Data types are flexible, including discrete and categorical data, such as soil or land-use maps, or continuous data as in digital elevation models, precipitation gradient maps, or pollutant concentration maps, etc.
Generally, two types of raster use in GIS and remote sensing application: a single band, or layer measure of a single characteristic of spatial area and multiple bands raster contains multiple spatially corresponding matrices of cell values representing the same spatial area. An example of a single-band raster data set is a digital elevation model (DEM). Each cell in a DEM contains only one value representing elevation of earth surface. Most satellite imagery has multiple bands, typically containing values within band of the electromagnetic spectrum.
In this exercise, we will use SRTM 90 Digital Elevation Model of Onondaga County, New York State which was downloaded from CGIAR-CSI.
Open With Rasterio¶
import rasterio as rio
from rasterio.plot import show
import pyproj
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
dem_path = "Data_02/Onondaga_DEM.tif"
dem=rio.open(dem_path)
#number of raster columns
dem.width
#number of raster bands
dem.count
#number of raster rows
dem.height
# All Metadata for the whole raster dataset
dem.meta
show(dem,cmap='RdYlGn')
from rasterio.plot import show_hist
show_hist(dem, bins=20, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram")
Open with GDAL¶
You can also open a raster layer with gdal
from osgeo import gdal
raster = gdal.Open(dem_path)
# Check type of the variable 'raster'
type(raster)
# Projection
raster.GetProjection()
data = raster.GetRasterBand(1).ReadAsArray()
plt.imshow(data,cmap='gray')
plt.show()
Write raster data¶
With Rasterio¶
We can export a raster file in python using the rasterio write() function.
Use rio.open() to create a new blank raster ‘template’. Then write the DEM numpy array to to that template using dst.write().
outTIF= "Data_02/Onondaga_DEM_NEW.tif"
dem_meta=dem.profile
with rio.open(outTIF, 'w', **dem_meta) as dst:
dst.write(dem.read(1), 1)
src = rio.open("Data_02/Onondaga_DEM_NEW.tif")
show(src)
With GDAL¶
Now we will how to write raster with GDAL driver. We will write a method to save the data in a format provided by the user. The procedure is fairly straightforward: we get a handler to a driver (e.g. a GeoTIFF or Erdas Imagine format), we create the output file (giving a filename, number of rows, columns, bands, the data type), and then add the relevant metadata (projection, geotransform, ...). We then select a band from the output and copy the array that we want to write to that band.
import gdal #
g = gdal.Open ("Data_02/Onondaga_DEM_NEW.tif" )
g.GetProjection()
# Get the x, y and number of bands from the original file
x_size, y_size, n_bands = g.RasterXSize, g.RasterYSize, g.RasterCount
data = g.ReadAsArray ()
# Get a handler to a driver
driver = gdal.GetDriverByName ( "HFA" )
# Next line creates the output dataset with
# 1. The filename ("Onondaga_DEM.img")
# 2. The raster size (x_size, y_size)
# 3. The number of bands
# 4.The data type (in this case, Byte.
# Other typical values might be gdal.GDT_Int16
# or gdal.GDT_Float32)
dataset_out = driver.Create( "Data_02/Onondaga_DEM_gdal.img", x_size, y_size, n_bands, \
gdal.GDT_Byte )
# Set the output geotransform by reading the input one
dataset_out.SetGeoTransform ( g.GetGeoTransform() )
# Set the output projection by reading the input one
dataset_out.SetProjection ( g.GetProjectionRef() )
# Now, get band # 1, and write our data array.
# Note that the data array needs to have the same type
# as the one specified for dataset_out
dataset_out.GetRasterBand ( 1 ).WriteArray ( data )
# This bit forces GDAL to close the file and write to it
dataset_out = None
!gdalinfo "Data_02/Onondaga_DEM_gdal.img"
Multi Band Raster¶
For loading multi-bands image in Python, we will use a sub-set of Landsat 8 multispectral image covering Onondaga county of New York state. This image was downloaded from USGS Earth Explore.
b2="Data_02/LC08_20150706_subset/LC08_L1TP_015030_20150716_20170226_01_T1_B2.TIF"
b3="Data_02/LC08_20150706_subset/LC08_L1TP_015030_20150716_20170226_01_T1_B3.TIF"
b4="Data_02/LC08_20150706_subset/LC08_L1TP_015030_20150716_20170226_01_T1_B4.TIF"
b5="Data_02/LC08_20150706_subset/LC08_L1TP_015030_20150716_20170226_01_T1_B5.TIF"
b2_r=rio.open(b2)
b3_r=rio.open(b3)
b4_r=rio.open(b4)
b5_r=rio.open(b5)
# Initialize subplots
import matplotlib.pyplot as plt
%matplotlib inline
fig, (ax1, ax2, ax3, ax4) = plt.subplots(ncols=4, nrows=1, figsize=(10, 4), sharey=True)
# Plot 4 bands
show(b2_r, ax=ax1)
show(b3_r, ax=ax2)
show(b4_r, ax=ax3)
show(b5_r, ax=ax4)
# Add titles
ax1.set_title("Band-2")
ax2.set_title("Band-3")
ax3.set_title("Band-4")
ax4.set_title("Band-5")
Stack Multi Band Imagery¶
Some remote sensing datasets are stored with each band in a separate file. However, often you want to use all of the bands together in your analysis. For example you need all of the bands together in the same file or “stack” in order to plot a color RGB image. EarthPy has a stack() function that allows you to take a set of .tif files that are all in the same spatial extent, CRS and resolution and either export them together a single stacked .tif file or work with them in Python directly as a stacked numpy array.
To begin using the EarthPy stack() function, import the needed packages and create an array to be plotted. Below you plot the data as continuous with a colorbar using the plot_bands() function.
import os
import numpy as np
# File manipulation
from glob import glob
import matplotlib.pyplot as plt
import geopandas as gpd
import rasterio as rio
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
You will use the glob function see the files in landsat data directory
glob ("Data_02/LC08_20150706_subset/*")
Now we will create a list for band (2,3,4 and 5)
multi_bands=glob( "Data_02/LC08_20150706_subset/*B[2:3:4:5]*.tif")
multi_bands
Now you have a list of all of the landsat bands in your landsat collections folder. You could chose to open each file individually using the rio.open (rasterio library) function.
Remember that Python uses 0 based indexing so band 3 is actually at index [2] not [3].
# Open a single band (NIR) and plot
with rio.open(multi_bands[3]) as src:
landsat_band4 = src.read()
ep.plot_bands(landsat_band4[0],
title="Landsat NIR (Band-5)",
scale=False)
plt.show()
Crerate Raster stack¶
landsat_multi_path = "Data_02/LC08_20150706_subset/landsat_multi.tif"
# This will create a new stacked raster with all bands save as landsat_multi.tif
land_stack, land_meta = es.stack(multi_bands,
landsat_multi_path)
Once we have stacked 4 bands, we can import it and work with it as we need to!
with rio.open(landsat_multi_path) as src:
landsat_multi = src.read()
# Plot all bands using earthpy
band_titles = ["Blue", "Green", "Red", "NIR"]
ep.plot_bands(landsat_multi,
title=band_titles, cbar=False)
plt.show()
Sources:¶
https://geohackweek.github.io/raster/04-workingwithrasters/ https://www.earthdatascience.org/workshops/gis-open-source-python/subtract-rasters-in-open-source-python/ https://www.e-education.psu.edu/geog489/node/2215 https://www.earthdatascience.org/courses/earth-analytics-python/lidar-raster-data/customize-matplotlib-raster-maps/ https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html https://automating-gis-processes.github.io/2016/Lesson7-read-raster.html https://www.neonscience.org/plot-neon-rgb-py https://automating-gis-processes.github.io/CSC/notebooks/L5/plotting-raster.html https://www.earthdatascience.org/courses/earth-analytics-python/multispectral-remote-sensing-in-python/naip-imagery-raster-stacks-in-python/