Since it takes quite a bit of work to use the open source geospatial libraries, we have started investigating creating a Python library that wraps the open source libraries into a more easy to use package.
The goals of this effort are to see if we can create a package that:
SpPy relies upon a number of open source libraries so you'll also need to make sure the following libraries are installed in your selected version of Python. We have had the most success using the Unofficial Windows Binaries for Python Extension Packages to install the libraries. For Python 3.x, the pip utility also appears to be working.
SpPy is available from GitHub. We recommend downloading the entire project folder and place it in a folder that is accessible to your Python scripts. Then, add a path to the SpPy folder to your Python scripts. The following lines of code will then import the key SpPy libraries and import the shapely geometry types:
# Open source spatial libraries
import sys
sys.path.append("C:/Temp/SPPy") # add the path to the SpPy libraries
# Open source spatial libraries
from shapely.geometry import Point, LineString, Polygon,MultiPolygon # classes to work with vector data
import shapely
# SP Libraries
import SPPlot
import SPVectors
import SPView
import SPRasters
import SPReferencing
The libraries have the following purposes:
Below we have setup some file paths that will be used in the examples below.
#################################### OutputFolderPath="C:\\Temp\\" CountriesFilePath="C:\\Temp\\Examples\\Data\\ne_110m_admin_0_countries.shp" TestRasterFilePath="C:\\Temp\\Examples\\Data\\test_raster.tif" ####################################
To get information from a specific data set, you'll want to create an SPDatasetVector object and then load a file into it. This allows you to get information on the data set, modify it's contents, and save it back out again. Below is an example where we have setup a Temp folder with the SpPy library including the Examples.
TheDataset=SPVectors.SPDatasetVector() #create a new layer
TheDataset.Load(CountriesFilePath) # load the contents of the layer
print("Type: "+format(TheDataset.GetType())) # get the type of data in the dataset
print("CRS: "+format(TheDataset.GetCRS())) # return the corodinate reference system
print("NumFeatures: "+format(TheDataset.GetNumFeatures())) # get the number of features in the dataset
print("Bounds: "+format(TheDataset.GetBounds())) # get the spatial bounds of the features
NumAttributes=TheDataset.GetNumAttributes()
print("NumAttributes: "+format(NumAttributes)) # return the definitions for the attributes (this function will change in the future)
Index=0
while (Index < NumAttributes):
print("Attribute: "+format(TheDataset.GetAttributeName(Index))+", Type: "+format(TheDataset.GetAttributeType(Index))+", Width: "+format(TheDataset.GetAttributeWidth(Index)))
Index+=1
# Save the result
TheDataset.Save(OutputFolderPath+"CountryClone.shp")
The following operations will take one or more data sets as the input and then typically return a transformed data set. The data set can be a path to a file, an existing data set object in memory, or a shapely geometry object. The advantage of these functions over other approaches is that you can have a single line perform an entire task. Then, you can use the result in the next function. This allows a series of geospatial transformations to be put into a single file in a very concise and yet readable way.
You can also perform "single-layer operations" on one layer at a time in just a few lines of code as below. For these functions, the "InputFile" can be a file path to a vector data set (typically a shapefile), an SPDatasetVector object, or a shapely geometry object. The example below shows that the operations can be "daisy-chained" without having to save files unless needed.
BufferedDataset=SPVectors.Buffer(CountriesFilePath,10) BufferedDataset.Save(OutputFolderPath+"Buffered_10.shp") CentroidDataset=SPVectors.Centroid(CountriesFilePath) CentroidDataset.Save(OutputFolderPath+"Centroid.shp") ConvexHullDataset=SPVectors.ConvexHull(CountriesFilePath) ConvexHullDataset.Save(OutputFolderPath+"ConvexHull.shp") SimplifiedDataset=SPVectors.Simplify(CountriesFilePath,10) SimplifiedDataset.Save(OutputFolderPath+"Simplify.shp")
The following functions work with loaded vector datasets and modify the attributes associated with each feature.
TheDataset=SPVectors.SPDatasetVector() #create a new layer
TheDataset.Load(CountriesFilePath) # load the contents of the layer
TheValue=TheDataset.GetAttributeValue("LABELRANK",0) # get the value of the atteribute "LABELRANK" for the feature in row 0
print("The Attribute="+format(TheValue))
TheDataset.SetAttributeValue("LABELRANK",0,0) # set the attribute "sum" for the feature in row 0
TheDataset.DeleteAttribute("scalerank") # Remove an attribute
TheDataset.AddAttribute("Test","str",254,"hi") # add a new attribute with the name, type, width, and default value
TheDataset.Save(OutputFolderPath+"AddAndRemoveAttributes.shp") # save the output
You can create selections of features in a dataset using boolean comparison function calls. The functions take the name of an attribute column and a value to compare each of the attribute values to. These functionswill return an array that contains True where the comparison is True for a feature's specified attribute and False otherwise. You can then use this array to subset the dataset to just the features that had True returned. Note that you can also perform operations on the selection arrays such as AND and OR to create more complex selection critieria.
TheDataset=SPVectors.SPDatasetVector() #create a new layer
TheDataset.Load(CountriesFilePath) # load the contents of the layer
Selection=TheDataset.SelectGreater("LABELRANK",2) # select the features with label rank greater than 2
# Below are four examples of selecting features
Selection=TheDataset.SelectGreaterThanOrEqual("LABELRANK",2) # select the features with label rank greater than or equal to 2
Selection=TheDataset.SelectLess("LABELRANK",2) # select the features with label rank less than or equal to 2.
Selection=TheDataset.SelectLessThanOrEqual("LABELRANK",2) # select the features with label rank less than 2.
Selection=TheDataset.SelectEqual("ADMIN","United States of America") # find the features that have the matching name
# Then, we can subset the dataset using a selection and write it out to a new shapefile
TheDataset.SubsetBySelection(Selection) # sub set the dataset to just contain features for the United States of America
TheDataset.Save(OutputFolderPath+"Selection.shp") # save the output
You can delete individual features using the index (row in the attribute table) of the feature.
TheDataset=SPVectors.SPDatasetVector() #create a new layer TheDataset.Load(CountriesFilePath) # load the contents of the layer # delete the first feature in the dataset TheDataset.DeleteFeature(0) # remove a feature based on its row index from 0 TheDataset.Save(OutputFolderPath+"OneFeatureDeleted.shp") # save the output
One of the features of SpPy is that you can work with geometrics, datasets in memory, or files. The code below creates a shapely geometry, then creates a new dataset and adds the geometry to the dataset, and then saves the dataset to a shapefile. We'll use these objects in the examples below. Note that when a feature is added, the attributes for the feature are set to 0 or blank values.
# Create a rectangular polygon Top=45 Bottom=18 Left=-125 Right=-100 BigBox=shapely.geometry.Polygon([(Left,Top), (Right,Top), (Right,Bottom), (Left,Bottom),(Left,Top)]) # Create a new layer (dataset) BigBoxDataset=SPVectors.SPDatasetVector() #create a new layer # add the rectangular polygon BigBoxDataset.AddFeature(BigBox) BigBoxDataset.Save(OutputFolderPath+"BigBox.shp") # save the output
The code below shows that you can use geometries, like BigBox created above, to perform overlay operations on spatial data.
# Crop a shapefile with a polygon NewLayer=SPVectors.Intersection(OutputFolderPath+"Selection.shp",BigBox) # perform an intersection on a geometry to get the new layer NewLayer.Save(OutputFolderPath+"IntersectionWithGeometry.shp") # Union the shapefile with the polygon NewLayer=SPVectors.Union(OutputFolderPath+"Selection.shp",BigBox) # perform a union on a geometry to get the new layer NewLayer.Save(OutputFolderPath+"UnionWithGeometry.shp") # save the output # Find the difference between the shapefile and the polygon (erase the area of the BigBox) NewLayer=SPVectors.Difference(OutputFolderPath+"Selection.shp",BigBox) # NewLayer.Save(OutputFolderPath+"DifferenceWithGeometry.shp") # save the output
You can also perform overlay operations on datasets.
NewLayer=SPVectors.Difference(OutputFolderPath+"USA.shp",BigBoxDataset) # NewLayer.Save(OutputFolderPath+"DifferenceWithDataset.shp") # save the output
Finally, you can perform operations on shapefiles as well.
NewLayer=SPVectors.Difference(OutputFolderPath+"USA.shp",OutputFolderPath+"BigBox.shp") NewLayer.Save(OutputFolderPath+"DifferenceWithShapefile.shp") # save the output
The code below will load a raster into memory and then return basic data for the raster. The Start of the file is the same except we need access to one of
TheDataset=SPRasters.SPDatasetRaster()
TheDataset.Load(TestRasterFilePath)
print("Width in pixels: "+format(TheDataset.GetWidthInPixels()))
print("Height in pixels: "+format(TheDataset.GetHeightInPixels()))
print("Num Bands: "+format(TheDataset.GetNumBands()))
print("Projection: "+format(TheDataset.GetProjection()))
print("Resolution (x,y): "+format(TheDataset.GetResolution()))
TheBounds=TheDataset.GetBounds()
print("TheBounds (XMin,YMin,YMax,XMax): "+format(TheBounds))
TheDataset.Save(OutputFolderPath+"CopiedRaster.tif")
The code below will create a raster of a given width and height.
WidthInPixels=100 HeightInPixels=100 TheDataset=SPRasters.SPDatasetRaster() TheDataset.SetWidthInPixels(WidthInPixels) TheDataset.SetHeightInPixels(HeightInPixels) TheDataset.SetType(gdal.GDT_Float32) TheBands=TheDataset.AllocateArray() TheBands=TheBands[0] Row=0 while (Row<HeightInPixels): Column=0 while (Column<WidthInPixels): TheBands[Row][Column]=Column Column+=1 Row+=1 print(TheBands) TheDataset.SetUTM(10) TheDataset.SetNorthWestCorner(400000,4000000) TheDataset.SetResolution(30,30) TheDataset.Save(OutputFolderPath+"NewRaster.tif")
The raster functions are very flexible and can take file paths to a raster, a raster object, or a scalar value. Mathematical, boolean, and comparison functions are all available in SpPy. A few examples are below and the rest are listed in the reference at the end of this page.
# Find pixels that are less than 10 using a shapefile and a function NewRaster=SPRasters.LessThan(TestRasterFilePath,10) # results in 1 for pixels in Raster1 are less than pixels in Raster2 # Find pixels that are less than 10 using a dataset and an overloaded operator TheRaster=SPRasters.SPDatasetRaster() # create the dataset TheRaster.Load(TestRasterFilePath) # load the data from a file NewRaster=TheRaster < 10 # Use the less than operator instead of a function NewRaster.Save(OutputFolderPath+"LessThan.tif") # save the resut to a file
To change the CRS for a dataset, we use the Proj4 parameters. Note that the projection engine for rasteres, GDAL, requires the raster to be written to a file so you need to specify the output file path when projecting rasters and this is a different function call.
# here is how you would specify the projection for Alber's Equal Area
ProjectionCode="aea"
Parameters={
"lon_0":0,
"lat_1":60,
"lat_2":40
}
# here is an example projecting the test raster back to geographic
ProjectionCode="longlat"
Parameters={
}
TheRaster=SPRasters.SPDatasetRaster()
TheRaster.Load(TestRasterFilePath)
NewLayer=SPReferencing.ProjectRaster(TheRaster,ProjectionCode,Parameters,OutputFolderPath+"Projected.tif")
The SPView class allows dataset to be rendred into a view which can then be saved to a raster file or displayed in a window. Since the rendering requires symobology, the dataset is placed inside a layer object first and then sent to the view.
Note that this portion of SpPy is in development (as most of it is ;-)
TheDataset=SPVectors.SPDatasetVector() #create a new layer TheDataset.Load(CountriesFilePath) # load the contents of the layer TheBounds=TheDataset.GetBounds() MinX=TheBounds[0] MinY=TheBounds[1] MaxX=TheBounds[2] MaxY=TheBounds[3] TheView=SPView.SPView(500,500) TheView.SetBounds(MinX,MaxX,MinY,MaxY) # TheLayer=SPVectors.SPLayerVector() TheLayer.SetDataset(TheDataset) TheLayer.Render(TheView) TheView.Save(OutputFolderPath+"RenderedVectors.png")
Below is an initial reference for the available functions.
The functions below are available for an SPDatasetVector object.
| Name | Description |
|---|---|
| Data Management | |
| Load(FilePath) | Load data from the specified file path |
| Save(FilePath) | Save the data set to the specified path |
| Type=GetType() | Returns the type of vector data from the OGC standards (e.g. "Polygon", "Point", "MultiPolygon", etc). |
| SetType(NewType) | Sets the type of vector data from the OGC standards (e.g. "Polygon", "Point", "MultiPolygon", etc). Should only be called before data has been added to the data set. |
| CRS=GetCRS() | Returns the coordinate reference system information (see https://fiona.readthedocs.io/en/latest/manual.html) |
| SetCRS(CRS) | |
| Attribute management | |
| Value=GetAttributeValue(AttributeName,Row) | Return the value of the attribute with the specific name and row index. |
| SetAttributeValue(AttributeName,Row,NewValue) | Set the value of the attribute with the specified name and row index to the NewValue. |
| Num=GetNumAttributes() | Returns the number of attribute columns defined for the data set. |
| Name=GetAttributeName(ColumnIndex) | Returns the attribute name for the specified attribute column |
| Type=GetAttributeType(ColumnIndex) | Returns the attribute type ("str","int""float") for the specified attribute column |
| Width=GetAttributeWidth(ColumnIndex) | Returns the width of the attribute (really its precision) |
| AddAttribute(Name,Type,Width=None,Default=None) | Add a new attribute column with the specified parameters |
| DeleteAttribute(Name) | Delete the attribute row of the specified name |
| RowIndex=GetAttributeColumn(Name) | Return an index to the attribute column that has the specified name |
| Selection=SelectEqual(AttributeName,Value) | Return a selection list with True where the attribute values are equal to the specified value and False otherwise |
| Selection=SelectGreater(AttributeName,Value) | Return a selection list with True where the attribute values are greater than to the specified value and False otherwise |
| Selection=SelectGreaterThanOrEqual(AttributeName,Value) | Return a selection list with True where the attribute values are greater than or equal to the specified value and False otherwise |
| Selection=SelectLess(AttributeName,Value) | Return a selection list with True where the attribute values are less than to the specified value and False otherwise |
| Selection=SelectLessThanOrEqual(AttributeName,Value) | Return a selection list with True where the attribute values are less than or equal to the specified value and False otherwise |
| SubsetBySelection(Selection) | Remove any features from the data set that have False in the selection array. |
| Feature/Geometry Management | These functions are provided to make it easy to obtain information from the features geometries |
| NumFeatures=GetNumFeatures() | Return the number of features/rows in the dataset |
| DeleteFeature(Row) | Remove the specified feature, including its attributes, from the data set. |
| AddFeature(TheGeometry,TheAttributes=None): | Add a new geometry to the data set. An optional list of attribute values can also be specified. |
| TheGeometry=GetGeometry(FeatureIndex) | Return the shapely geometry object associated with the feature, if any |
| SplitFeatures() | Breaks up features that are collections or contain muiltple geometries (i.e. "multi") into individual features. |
| Bounds=GetBounds() | Return the spatial bounds of all the features. |
| Area=GetFeatureArea(FeatureIndex) | Returns the area of the feature |
| Bounds=GetFeatureBounds(FeatureIndex) | Returns the spatial bounds of the feature |
| Length=GetFeatureLength(FeatureIndex) | Returns the length of all the line segments in the feature |
| EmptyFlag=IsFeatureEmpty(FeatureIndex) | Returns True if the feature does not contain any coordinates |
| ValidFlag=IsFeatureValid(FeatureIndex) | Returns True if the feature contains a valid geometry (e.g. no overlapping segments) |
The following functions are available in the SPVector library. These functions take one or more inputs that can be paths to files, data sets loaded into memory, or shapely geometries.
| Single-Dataset Operations | |
| OutputDataset=Buffer(InputDataset,Amount) | |
| OutputDataset=Simplify(InputDataset,Tolerance, PreserveTopology=True) | |
| OutputDataset=ConvexHull(InputDataset) | |
| OutputDataset=Centroid(InputDataset) | |
| Overlay Operations | |
| OutputDataset=Intersection(InputDataset1,InputDataset2) | |
| OutputDataset=Union(InputDataset1,InputDataset2) | |
| OutputDataset=Difference(InputDataset1,InputDataset2) | |
| OutputDataset=SymmetricDifference(InputDataset1,InputDataset2) |
Provides functions to work with raster data sets. The operations that work on rasters can take a path to a file, an existing raster object, or a scalar value for operators.
| Name | Description |
|---|---|
| RasterDataset operations | |
| RasterDataset.Load(FilePath) | Loads the data from a file into the dataset |
| RasterDataset.Save(FilePath) | Saves the data from the dataset into a file |
| WidthInPixels=RasterDataset.GetWidthInPixels() | Retursn the number of columns of pixels in the data |
| RasterDataset.SetWidthInPixels(NewWidth) | Sets the overall width of the raster in pixels. Only should be called if the raster data has not been allocated yet. |
| HeightInPixels=RasterDataset.GetHeightInPixels() | Returns the overall width of the raster in pixels. Only should be called if the raster data has not been allocated yet. |
| RasterDataset.SetHeightInPixels(NewHeight) | Sets the number of rows of pixels in the data |
| NumBands=RasterDataset.GetNumBands() | Returns the number of bands of data in the raster (i.e. gray=1, RGB=3). |
| CRS=RasterDataset.GetProjection() | Note: Needs to change to GetCRS() |
| X,Y=RasterDataset.GetResolution() | Returns the resolution (spatial reference units per pixels) as a tuple |
| RasterDataset.SetResolution((X,Y)) | |
| Bounds=RasterDataset.GetBounds() | |
| RasterDataset.SetNorthWestCorner(X,Y) | |
| Band=RasterDataset.GetBand(BandIndex) | Returns an array with the band data for the specified index (0 for first, 1 for second, etc.) |
| Bands=RasterDataset.GetBands() | Return an array of arrays with the data for all the bands. |
| RasterDataset.SetBands(Bands) | Sets an array of arrays with the data for all the bands. This function should be used carefully as the dimensions and type of the band data must match the settings in the object. |
| Min,Max=RasterDataset.GetMinMax(BandIndex) | REturns the minimum and maximum values for the data in the specific bands. |
| GDALType=RasterDataset.GetType() | Returns the type of sample values based on GDAL definitions: gdal.GDT_Byte, gdal.GDT_UInt16,gdal.GDT_Int16,gdal.GDT_UInt32,gdal.GDT_Int32,gdal.GDT_Float32,dal.GDT_Float64,gdal.GDT_CFloat64 Note: Do we support all these types? |
| RasterDataset.SetType(New) | Sets the type of sample values based on GDAL definitions (see above) |
| RasterDataset.AllocateArray() | Allocates the pixel data for the current data set |
SPRaster Public Functions
The following functions are available in the SPRaster library as single-line operations. The inputs can be paths to files or data sets loaded into memory.
| Math Operations | |
OutputRaster=Add(Input1,Input2) Outputraster=Input1+Input2 |
Adds the pixels in the two rasters together or add a scalar value to the pixels in the first raster. |
OutputRaster=Subtract(Input1,Input2) OutputRaster=Input1 - Input2 |
Subtract the pixels in the two rasters from one another or subtract a scalar value from the raster. |
OutputRaster=Multiply(Input1,Input2) OutputRaster=Input1 * Input2 |
Mutiply the pixels in the two rasters together or mutiply a scalar value to the pixels in the raster. |
OutputRaster=Divide(Input1,Input2) OutputRaster=Input1 / Input2 |
Divide the raster Input1 by the raster or scalar value in Input2. |
| Comparisons | |
OutputRaster=Equal(Input1,Input2) OutputRaster=Input1==Input2 |
Creates a new raster with True (1) for pixels that are equal in both rasters and 0 otherwise. |
OutputRaster=LessThan(Input1,Input2) OutputRaster=Input1<Input2 |
|
OutputRaster=GreaterThan(Input1,Input2) OutputRaster=Input1>Input2 |
|
OutputRaster=LessThanOrEqual(Input1,Input2) OutputRaster=Input1<=Input2 |
|
OutputRaster=GreaterThanOrEqual(Input1,Input2) OutputRaster=Input1>=Input2 |
|
| OutputRaster=Minimum(Input1,Input2) | Returns a new raster where each pixel is the minimum value of the two overlapping pixels in the inputs. |
| OutputRaster=Maximum(Input1,Input2) | |
| Boolean Operations | |
OutputRaster=And(Input1,Input2) |
Returns a new raster where each pixel is the logical result of anding the overlapping pixels in the inputs. |
| OutputRaster=Or(Input1,Input2) | |
| OutputRaster=Not(Input1) | |
| Conversions | |
| OutputRaster=RoundInteger(Input1) | Returns a new raster where the pixels have been rounded to the nearest integer. |
| OutputRaster=RoundFix(Input1) | Round to nearest integer towards zero. |
| OutputRaster=RoundFloor(Input1) | Returns a new raster where the pixels have been truncated to the lowest value (i.e. 0.5 -> 0, 0.6 -> 0) |
| OutputRaster=RoundCeiling(Input1) | Returns a new raster where the pixels have been truncated to the highest nvalue (i.e. 0.5 -> 1, 0.6 -> 1) |
| OutputRaster=RoundTruncate(Input1) | Returns a new raster where the pixels have been truncated (i.e. 0.5 -> 0, 0.6 -> 0) |
| Exponential functions | |
| OutputRaster=NaturalLog(Input1) | |
| OutputRaster=Log(Input1) | |
| OutputRaster=Exponent(Input1) | Returns a new raster with pixels that are the result of raising the value e to the value in the input raster. |
| OutputRaster=Power(Input1,Power) | Returns a new raster with pixels that are the result of raising the values in the input by the specified Power (i.e. New=Old ^Power). |
| OutputRaster=Square(Input1) | |
| OutputRaster=SquareRoot(Input1) | |
| OutputRaster=AbsoluteValue(Input1) | |
| Conversions | |
| OutputVector=Polygonize(Input1) | Need to finish this jjg |
This class provides resampling for raster data
| Name | Description |
|---|---|
| NewDataset=Zoom(Input1,ZoomFactor) | |
| ExtractByPixels(InputRasterDataset,StartColumn,StartRow,EndColumn,EndRow) |
Provides functions to project datasets.
| Name | Description |
|---|---|
| NewDataset=SPReferencing.Project(InputDataset,ProjectionCode,Parameters,LatNewPole=90) | Projecs the dataset to the desired spatial reference |
© Copyright 2018 HSU - All rights reserved.