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:
The key to this effort is to see if we can create a design for the package that meets the goals. This effort is underway and a prototype of the code is available for you to evaluate. Your feedback on the design of the approach is critical at this stage. If successful, we'll look at flushing out the package with additional code and setting up formal testing and documentation.
To get started, download the SpPy (spatial pie) package here 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 from shapely.geometry import Point, LineString, Polygon,MultiPolygon # SP Libraries import SPPlot import SPVectors import SPView import SPRasters import SPReferencing
The libraries have the following purposes:
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.
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(InputDataset,10) CentroidDataset=SPVectors.Centroid(BufferedDataset) ConvexHullDataset=SPVectors.ConvexHull(CentroidDataset) SimplifiedDataset=SPVectors.Simplify(ConvexHullDataset,10) SimplifiedDataset.Save(OutputFolderPath+"Simplify.shp")
Finally, there are overlay operations that can be executed in just a couple lines of code.
# Create the bounding polygon Top=45 Bottom=20 Left=-125 Right=-40 BoundingPoly=shapely.geometry.Polygon([(Left,Top), (Right,Top), (Right,Bottom), (Left,Bottom),(Left,Top)]) # Crop a shapefile with a polygon NewLayer=SPVectors.Intersection(InputFile,BoundingPoly) # perform an intersection on a geometry to get the new layer NewLayer.Save(OutputFolderPath+"Intersection.shp") # Union the shapefile with the polygon NewLayer=SPVectors.Union(InputFile,BoundingPoly) # perform a union on a geometry to get the new layer NewLayer.Save(OutputFolderPath+"Union.shp") # save the output # Find the difference between the shapefile and the polygon NewLayer=SPVectors.Difference(InputFile,BoundingPoly) # NewLayer.Save(OutputFolderPath+"Difference.shp") # save the output
You can easily add and remove features using the commands below. When a feature is added, the attributes for the feature are set to 0 or blank values.
# add a square geometry in at 0,0 TheGeometry=shapely.geometry.Polygon([(-10,10), (10,10), (10,-10), (-10,-10),(-10,10)]) TheDataset.AddFeature(TheGeometry) # delete the first feature in the dataset TheDataset.DeleteFeature(0) # remove a feature based on its row index from 0
The following functions work with loaded vector datasets and modify the attributes associated with each feature.
TheDataset.GetAttributeValue("LABELRANK",0) # get the value of the atteribute "LABELRANK" for the feature in row 0
TheDataset.SetAttributeValue("sum",0,"test2") # set the attribute "sum" for the feature in row 0
TheDataset.DeleteAttribute("ADM0_DIF") # Remove an attribute
TheDataset.AddAttribute("testy","str",254,"hi") # add a new attribute with the name, type, width, and default value
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.
Selection=TheDataset.SelectGreater("LABELRANK",2) # select the features with label rank greater than 2
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
TheDataset.SubsetBySelection(Selection) # sub set the dataset to just contain features for the United States of America
The code below will load a raster into memory and then return basic data for the raster.
TheDataset =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))
#TheBand=TheDataset.GetBandAsArray(1)
#print("TheBands: "+format(TheBand))
TheDataset.Save(TempFolderPath+"CopiedRaster.tif")
WidthInPixels=100 HeightInPixels=100 TheDataset=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(TempFolderPath+"NewRaster.tif")
Result=SPRaster.LessThan(TestRasterFilePath,TestRasterFilePath2)
To change the CRS for a dataset, we use the Proj4 parameters.
ProjectionCode="aea"
Parameters={
"lon_0":0,
"lat_1":60,
"lat_2":40
}
NewLayer=SPReferencing.Project(InputDataset,ProjectionCode,Parameters)
NewLayer.Save(OutputFolderPath+"Projected.shp")
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.
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) # jjg - should be able to set the bounds directly 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) | Adds the pixels in the two rasters together. |
| OutputRaster=Subtract(Input1,Input2) | |
| OutputRaster=Multiply(Input1,Input2) | |
| OutputRaster=Divide(Input1,Input2) | |
| Comparisons | |
| OutputRaster=Equal(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=GreaterThan(Input1,Input2) | |
| OutputRaster=LessThanOrEqual(Input1,Input2) | |
| OutputRaster=GreaterThanOrEqual(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.