Read Field Values of a Shapefile Gdal
Shapefile (Vector) Handling¶
This section introduces geospatial analysis of shapefiles with gdal, ogr, and osr. For interactive reading and executing code blocks 
and discover geo01-shp.ipynb, or install Python and JupyterLab locally.
Use flusstools
The core functions featured in this department are besides implemented in flusstools. To use those functions, brand sure flusstools is installed and import information technology as follows: from flusstools import geotools . Some of the functions shown in this tutorial tin then exist used with geotools.function_name() .
Load an Existing Shapefile¶
OSGeo's ogr module is an splendid source for handling shapefiles. Later importing the libraries, the correct driver for treatment shapefiles ( "ESRI Shapefile" ) with Python needs to exist chosen. Next, with the commuter object ( ogr.GetDriverByName("SHAPEFILE") ), we can open up (instantiate) a shapefile (object with shp_driver.Open("SHAPEFILE") ), which contains layer information. It is precisely the layer data (i.e., references to shapefile attributes) that we desire to work with. Therefore we take to instantiate a shapefile layer object with shp_dataset.GetLayer() .
from osgeo import ogr shp_driver = ogr . GetDriverByName ( "ESRI Shapefile" ) shp_dataset = shp_driver . Open ( "geodata/shapefiles/cities.shp" ) shp_layer = shp_dataset . GetLayer ()
Practice not use from gdal import ogr
Importing ogr with from gdal import ogr is deprecated since GDAL v3. Therefore, ogr must be imported from osgeo .
Create a New Shapefile¶
ogr too enables creating a new bespeak, line, or polygon shapefile. The post-obit code block defines a function for creating a shapefile, where the optional keyword argument overwrite is used to command whether an existing shapefile with the same name should exist overwritten (default: True ). The command shp_driver.CreateDataSource(SHP-FILE-DIR) creates a new shapefile and the residue of the function adds a layer to the shapefile if the optional keyword arguments layer_name and layer_type are provided. Both optional keywords must be strings, where layer_name can be any proper name for the new layer. layer_type must be either "betoken" , "line" , or "polygon" to create a betoken, (poly)line, or polygon shapefile, respectively. The role uses the geometry_dict dictionary to assign the correct ogr.SHP-TYPE to the layer. There are more options for extending the create_shp(...) function listed on pcjerick's Github pages.
from osgeo import ogr import bone def create_shp ( shp_file_dir , overwrite = Truthful , * args , ** kwargs ): """ :param shp_file_dir: STR of the (relative shapefile directory (ends on ".shp") :param overwrite: [optional] BOOL - if Truthful, existing files are overwritten :kwarg layer_name: [optional] STR of the layer_name - if None: no layer will be created (max. 13 chars) :kwarg layer_type: [optional] STR ("indicate, "line", or "polygon") of the layer_name - if None: no layer will be created :output: returns an ogr shapefile layer """ shp_driver = ogr . GetDriverByName ( "ESRI Shapefile" ) # check if output file exists if aye delete it if bone . path . exists ( shp_file_dir ) and overwrite : shp_driver . DeleteDataSource ( shp_file_dir ) # create and return new shapefile object new_shp = shp_driver . CreateDataSource ( shp_file_dir ) # create layer if layer_name and layer_type are provided if kwargs . get ( "layer_name" ) and kwargs . become ( "layer_type" ): # create dictionary of ogr.SHP-TYPES geometry_dict = { "point" : ogr . wkbPoint , "line" : ogr . wkbMultiLineString , "polygon" : ogr . wkbMultiPolygon } # create layer try : new_shp . CreateLayer ( str ( kwargs . get ( "layer_name" )), geom_type = geometry_dict [ str ( kwargs . get ( "layer_type" ) . lower ())]) except KeyError : print ( "Error: Invalid layer_type provided (must be 'betoken', 'line', or 'polygon')." ) except TypeError : print ( "Error: layer_name and layer_type must be string." ) except AttributeError : print ( "Error: Cannot access layer - opened in other programme?" ) return new_shp
The create_shp function is as well provided with flusstools ( flusstools.geotools.create_shp() ) and aids to create a new shapefile (brand certain to get the directory right):
a_new_shp_file = create_shp ( r "" + bone . getcwd () + "/geodata/shapefiles/new_polygons.shp" , layer_name = "basemap" , layer_type = "polygon" ) # release data source a_new_shape_file = None
Important
A shapefile name may non accept more than thirteen characters and a field name may not take more than x characters (read more in Esri'southward shapefile docs).
Shapefiles can as well be created and fatigued in QGIS and the following figures guide through the procedure of creating a polygon shapefile. We will not demand the resulting shapefile in this department anymore, simply later for making information technology interact with raster datasets.
The first step to making a shapefile with QGIS is patently to run QGIS and create a new project. The post-obit example uses water depth and period velocity raster information every bit groundwork information to delineate the and so-called morphological unit of slackwater. Both the water depth and flow velocity rasters are office of the River Architect sample datasets (precisely located at RiverArchitect/SampleData/01_Conditions/2100_sample/ ). After downloading the sample data, they can be imported in QGIS by dragging the files from the Browser console into the Layers console. Then:
Get and Set Shapefile Projections¶
The terminology used in the .prj files of a shapefile corresponds to the definitions in the geospatial information section. In Python, information about the coordinate system is bachelor through shp_layer.GetSpatialRef() of the ogr library:
shp_srs = shp_layer . GetSpatialRef () print ( shp_srs )
GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, Dominance["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], Unit of measurement["degree",0.0174532925199433, Dominance["EPSG","9122"]], Axis["Latitude",NORTH], AXIS["Longitude",Due east], Say-so["EPSG","4326"]]
This GEOGCS definition of the above shapefile corresponds to Esri's well-known text (WKT). Since the shapefile format was developed by Esri, Esri's WKT (esriwkt) format must be used in .prj files. The Open up Geospatial Consortium (OGC) uses a unlike well-known text in their EPSG:XXXX definitions (e.1000., available at spatialreference.org).
GEOGCS [ "WGS 84" , DATUM [ "WGS_1984" , SPHEROID [ "WGS84" , 6378137 , 298.257223563 , AUTHORITY [ "EPSG" , "7030" ]] , AUTHORITY [ "EPSG" , "6326" ]] , PRIMEM [ "Greenwich" , 0 , Authorisation [ "EPSG" , "8901" ]] , Unit [ "degree" , 0.01745329251994328 , AUTHORITY [ "EPSG" , "9122" ]] , Potency [ "EPSG" , "4326" ]] To redefine or newly define the coordinate organization of a shapefile, we can use spatialreference.org through Python's default urllib library.
Notation
Running the following code block requires an net connection.
import urllib # role to get spatialreferences with epsg code def get_esriwkt ( epsg ): # usage get_epsg_code(4326) try : with urllib . asking . urlopen ( "http://spatialreference.org/ref/epsg/ {0} /esriwkt/" . format ( epsg )) every bit response : render str ( response . read ()) . strip ( "b" ) . strip ( "'" ) except Exception : pass try : with urllib . request . urlopen ( "http://spatialreference.org/ref/sr-org/epsg {0} -wgs84-web-mercator-auxiliary-sphere/esriwkt/" . format ( epsg )) equally response : return str ( response . read ()) . strip ( "b" ) . strip ( "'" ) # sr-org codes are available at "https://spatialreference.org/ref/sr-org/{0}/esriwkt/".format(epsg) # for case EPSG:3857 = SR-ORG:6864 -> https://spatialreference.org/ref/sr-org/6864/esriwkt/ = EPSG:3857 except Exception : impress ( "Mistake: Could not discover epsg lawmaking on spatialreference.org. Returning default WKT(epsg=4326)." ) return 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],Unit of measurement["Caste",0.017453292519943295],UNIT["Meter",1]]'
This function tin be used to create a new projection file:
# open the hypy-area shapefile shp_file = "hypy-area" # create new .prj file for the shapefile (.shp and .prj must take the same name) with open ( "geodata/shapefiles/ {0} .prj" . format ( shp_file ), "west" ) equally prj : epsg_code = get_esriwkt ( 4326 ) prj . write ( epsg_code ) impress ( "Wrote projection file : " + epsg_code )
Wrote projection file : GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],Unit of measurement["Degree",0.017453292519943295]]
An offline alternative for generating a .prj file is the osr library that comes along with gdal :
from osgeo import osr def get_wkt ( epsg , wkt_format = "esriwkt" ): default = 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Caste",0.017453292519943295],UNIT["Meter",1]]' spatial_ref = osr . SpatialReference () try : spatial_ref . ImportFromEPSG ( epsg ) except TypeError : print ( "ERROR: epsg must be integer. Returning default WKT(epsg=4326)." ) return default except Exception : print ( "Error: epsg number does not be. Returning default WKT(epsg=4326)." ) return default if wkt_format == "esriwkt" : spatial_ref . MorphToESRI () # return a nicely formatted WKT string (alternatives: ExportToPCI(), ExportToUSGS(), or ExportToXML()) return spatial_ref . ExportToPrettyWkt ()
Transform (Re-projection) a Shapefile¶
A re-projection may be needed, for instance, to utilise a shapefile in EPSG:4326 (e.g., created with QGIS) in a web-GIS application (due east.1000., open street maps) that typically uses EPSG:3857 . To utilize a different projection to geometric objects of a shapefile, it is not enough to simply rewrite the .prj file. Therefore, the following example shows the re-project of the countries.shp shapefile from the Natural Earth quick showtime kit. The following workflow performs the reprojection:
-
Save shapefile to transform is located in the subdirectory
geodata/shapefiles/countries.shpand open it in the Python script as described above. -
Read and identify the spatial reference system used in the input shapefile:
-
Create a spatial reference object with
in_sr = osr.SpatialReference(str(shapefile.GetSpatialRef())). -
Detect the spatial reference organization in EPSG format with
AutoIdentifyEPSG(). -
Assign the EPSG-formatted spatial reference organization to the spatial reference object of the input shapefile (
ImportFromEPSG(int(in_sr.GetAuthorityCode(None)))).
-
-
Create the output spatial reference with
out_sr = osr.SpatialReference()and use the target EPSG code (out_sr.ImportFromEPSG(3857)). -
Create a coordinate transformation object (
coord_trans = osr.CoordinateTransformation(in_sr, out_sr)) that enables re-projecting geometry objects subsequently. -
Create the output shapefile, which will correspond to a copy of the input shapefile (employ the above-defined
create_shpfunction withlayer_name="basemap"andlayer_type="line"). -
Copy the field names and types of the input shapefile:
-
Read the attribute layer from the input file's layer definitions with
in_lyr_def = in_shp_lyr.GetLayerDefn() -
Iterate through the field definitions and append them to the output shapefile layer (
out_shp_lyr)
-
-
Iterate through the geometry features in the input shapefile:
-
Employ the new (output) shapefile'south layer definitions (
out_shp_lyr_def = out_shp_lyr.GetLayerDefn()) to append transformed geometry objects later. -
Define an iteration variable
in_featureas an example ofin_shp_lyr.GetNextFeature. -
In a
whileloop, instantiate every geometry (geometry = in_feature.GetGeometryRef()) in the input shapefile, transform thegeometry(geometry.Transform(coord_trans)), convert it to anogr.Characteristic()with theSetGeometry(geometry)method, re-create field definitions (nestedfor-loop), and suspend the new feature to the output shapefile layer (out_shp_lyr_def.CreateFeature(out_feature)). -
At the end of the
while-loop, look for the side by side feature in the input shapefile'due south attributes within_feature = in_shp_lyr.GetNextFeature()
-
-
Unlock (release) all layers and shapefiles by overwriting the objects with
None(nothing is written to whatsoever file every bit long as these variables exist!). -
Assign the new projection EPSG:3857 using the in a higher place-defined
get_wktfunction.
from osgeo import ogr from osgeo import osr shp_driver = ogr . GetDriverByName ( "ESRI Shapefile" ) # open input shapefile and layer in_shp = shp_driver . Open ( r "" + os . path . abspath ( '' ) + "/geodata/shapefiles/countries.shp" ) in_shp_lyr = in_shp . GetLayer () # become input SpatialReference in_sr = osr . SpatialReference ( str ( in_shp_lyr . GetSpatialRef ())) # auto-detect epsg in_sr . AutoIdentifyEPSG () # assign input SpatialReference in_sr . ImportFromEPSG ( int ( in_sr . GetAuthorityCode ( None ))) # create SpatialReference for new shapefile out_sr = osr . SpatialReference () out_sr . ImportFromEPSG ( 3857 ) # create a CoordinateTransformation object coord_trans = osr . CoordinateTransformation ( in_sr , out_sr ) # create output shapefile and get layer out_shp = create_shp ( r "" + os . path . abspath ( '' ) + "/geodata/shapefiles/countries-spider web.shp" , layer_name = "basemap" , layer_type = "line" ) out_shp_lyr = out_shp . GetLayer () # await up layer (features) definitions in input shapefile in_lyr_def = in_shp_lyr . GetLayerDefn () # re-create field names of input layer attribute tabular array to output layer for i in range ( 0 , in_lyr_def . GetFieldCount ()): out_shp_lyr . CreateField ( in_lyr_def . GetFieldDefn ( i )) # instantiate feature definitions object for output layer (currently empty) out_shp_lyr_def = out_shp_lyr . GetLayerDefn () # iteratively suspend all input features in new project in_feature = in_shp_lyr . GetNextFeature () while in_feature : # go the input geometry geometry = in_feature . GetGeometryRef () # re-projection (transform) geometry to new system geometry . Transform ( coord_trans ) # create new output characteristic out_feature = ogr . Characteristic ( out_shp_lyr_def ) # assign in-geometry to output feature and copy field values out_feature . SetGeometry ( geometry ) for i in range ( 0 , out_shp_lyr_def . GetFieldCount ()): out_feature . SetField ( out_shp_lyr_def . GetFieldDefn ( i ) . GetNameRef (), in_feature . GetField ( i )) # add together the feature to the shapefile out_shp_lyr . CreateFeature ( out_feature ) # ready next iteration in_feature = in_shp_lyr . GetNextFeature () # release shapefiles and layers in_shp = None in_shp_lyr = None out_shp = None out_shp_lyr = None # create .prj file for output shapefile for web application references with open ( r "" + os . path . abspath ( '' ) + "/geodata/shapefiles/countries-spider web.prj" , "w+" ) as prj : prj . write ( get_wkt ( 3857 ))
Challenge
Re-write the above lawmaking cake in a re_project(shp_file, target_epsg) function.
In the example of success, the lawmaking sequence in_sr.AutoIdentifyEPSG() returns 0 (i.e., it successfully recognized the EPSG number), just unfortunately, many EPSG numbers are non known to AutoIdentifyEPSG() . In the case that AutoIdentifyEPSG() did non function propperly, the method does non render 0 , but 7 , for example. A workaround for the limited functionality of srs.AutoIdentifyEPSG() is srs.FindMatches . srs.FindMatches returns a matching srs_match from a larger database, which is somewhat nested; for instance, utilize:
matches = srs . FindMatches () Then, matches looks similar this: [(osgeo.osr.SpatialReference, INT)] . Therefore, a complete workaround for srs.AutoIdentifyEPSG() (or in_sr.AutoIdentifyEPSG() in the code block higher up) looks like this:
# set epsg and create spatial reference object epsg = 3857 srs = osr . SpatialReference () srs . ImportFromEPSG ( epsg ) # identify spatial reference auto_detect = srs . AutoIdentifyEPSG () if auto_detect is not 0 : srs = srs . FindMatches ()[ 0 ][ 0 ] # Find matches returns listing of tuple of SpatialReferences srs . AutoIdentifyEPSG () # Re-perform auto-identification
Add Fields and Point Features to a Shapefile¶
A shapefile feature can exist a bespeak, a line, or a polygon, which has field attributes (e.one thousand., "id"=one to describe that this is polygon number i or associated to an id block 1). Field attributes tin be more than but an IDentifier and include, for example, the polygon surface area or city labels as in the to a higher place-shown example illustrating cities.shp ( shp_driver.Open("geodata/shapefiles/cities.shp") ).
To create a signal shapefile, we tin can use the above-introduced create_shp function (or flusstools.geotools.create_shp(shp_file_dir="") ) and set the projection with the get_epsg_code() function (also higher up introduced). The following lawmaking cake shows the usage of both functions to create a river.shp indicate shapefile that contains three points located at three rivers in cardinal Europe. The code block also illustrates the creation of a field in the attribute table and the creation of three betoken features. Here is how the code works:
-
The shapefile is located in the
rivers_ptsvariable. Note that thelayer_typealready determines the type of geometries that can be used in the shapefile. For instance, adding a line or polygon feature to anogr.wkbPointlayer will upshot in anError 1message. -
The
basemap(layer) is attributed to the variablelyr = river_pts.GetLayer(). -
A string-type field is added and appended to the attribute table:
-
instantiate a new field with
field_gname = ogr.FieldDefn("FIELD-NAME", ogr.OFTString)(call back: the field name may not have more than 10 characters!) -
append the new field to the shapefile with
lyr.CreateField(field_gname) -
other field types than
OFTStringcan exist:OFTInteger,OFTReal,OFTDate,OFTTime,OFTDateTime,OFTBinary,OFTIntegerList,OFTRealList, orOFTStringList.
-
-
Add three points stored in
pt_names = {RIVER-NAME: (x-coordinate, y-coordinate)}in a loop over the dictionary keys:-
for every new point, create a characteristic every bit a kid of the layer defintions with
characteristic = ogr.Characteristic(lyr.GetLayerDefn()) -
set the value of the field proper noun for every indicate with
feature.SetField(FIELD-NAME, FIELD-VALUE) -
create a string of the new point in WKT format with
wkt = "POINT(10-COORDINATE Y-COORDINATE)" -
catechumen the WKT-formatted point into a point-blazon geometry with
signal = ogr.CreateGeometryFromWkt(wkt) -
set the new point every bit the new feature's geometry with
feature.SetGeometry(point) -
suspend the new feature to the layer with
lyr.CreateFeature(feature)
-
-
Unlock (release) the shapefile by overwriting the
lyrandriver_ptsvariable withNone.
Important
The operations are not written to the shapefile if the lyr and river_pts objects are not overwritten with None .
shp_dir = r "" + os . path . abspath ( '' ) + "/geodata/shapefiles/rivers.shp" river_pts = create_shp ( shp_dir , layer_name = "basemap" , layer_type = "point" ) # create .prj file for the shapefile for spider web application references with open up ( shp_dir . split up ( ".shp" )[ 0 ] + ".prj" , "w+" ) as prj : prj . write ( get_esriwkt ( 3857 )) # go basemap layer lyr = river_pts . GetLayer () # add string field "rivername" field_gname = ogr . FieldDefn ( "rivername" , ogr . OFTString ) lyr . CreateField ( field_gname ) # names and coordinates of central European union rivers in EPSG:3857 WG84 / Pseudo-Mercator pt_names = { "Aare" : ( 916136.03 , 6038687.72 ), "Ain" : ( 623554.12 , 5829154.69 ), "Inn" : ( 1494878.95 , 6183793.83 )} # add together the three rivers equally points to the basemap layer for north in pt_names . keys (): # create Characteristic as child of the layer pt_feature = ogr . Characteristic ( lyr . GetLayerDefn ()) # ascertain value north (river) in the rivername field pt_feature . SetField ( "rivername" , n ) # use WKT format to add a point geometry to the Feature wkt = "Betoken( %f %f )" % ( float ( pt_names [ n ][ 0 ]), float ( pt_names [ n ][ 1 ])) point = ogr . CreateGeometryFromWkt ( wkt ) pt_feature . SetGeometry ( bespeak ) # append the new characteristic to the basement layer lyr . CreateFeature ( pt_feature ) # release files lyr = None river_pts = None
The resulting rivers.shp shapefile tin be imported in QGIS forth with a DEM from the Natural Globe quick start kit.
Multiline (Polyline) Shapefile¶
Similar to the procedure for creating and calculation points to a new indicate shapefile, a (multi) line (or polyline) can be added to a shapefile. The create_shp function creates a multi-line shapefile when the layer type is defined as "line" ( flusstools.geotools.create_shp(shp_file_dir="", layer_type="line") ). The coordinate system is created with the higher up-defined get_gps_code function.
Tip
The term multi-line is used in OGC and ogr , while polyline is used in Esri GIS environments.
The following code block uses the coordinates of cities along the Rhine River, stored in a dictionary named station_names . The metropolis names are not used and only the coordinates are appended with line.AddPoint(X, Y) . As earlier, a field is created to give the river a name. The actual line feature is again created as a child of the layer with line_feature = ogr.Feature(lyr.GetLayerDefn()) . Running this code block produces a line that approximately follows the Rhine River between France and Germany.
shp_dir = r "" + os . path . abspath ( '' ) + "/geodata/shapefiles/rhine_proxy.shp" rhine_line = create_shp ( shp_dir , layer_name = "basemap" , layer_type = "line" ) # create .prj file for the shapefile for spider web application references with open ( shp_dir . split ( ".shp" )[ 0 ] + ".prj" , "west+" ) as prj : prj . write ( get_wkt ( 3857 )) # go basemap layer lyr = rhine_line . GetLayer () # coordinates for EPSG:3857 WG84 / Pseudo-Mercator station_names = { "Basel" : ( 844361.68 , 6035047.42 ), "Kembs" : ( 835724.27 , 6056449.76 ), "Breisach" : ( 842565.32 , 6111140.43 ), "Rhinau" : ( 857547.04 , 6158569.58 ), "Strasbourg" : ( 868439.31 , 6203189.68 )} # create line object and add points from station names line = ogr . Geometry ( ogr . wkbLineString ) for stn in station_names . values (): line . AddPoint ( stn [ 0 ], stn [ 1 ]) # create field named "rives" field_name = ogr . FieldDefn ( "river" , ogr . OFTString ) lyr . CreateField ( field_name ) # create feature, geometry, and field entry line_feature = ogr . Characteristic ( lyr . GetLayerDefn ()) line_feature . SetGeometry ( line ) line_feature . SetField ( "river" , "Rhine" ) # add together feature to layer lyr . CreateFeature ( line_feature ) lyr = None rhine_line = None
The resulting rhine_proxy.shp shapefile tin be imported in QGIS forth with a DEM and the cities signal shapefile from the Natural Earth quick first kit.
Polygon Shapefile¶
Polygons are surface patches that can exist created indicate-past-bespeak, line-by-line, or from a "Multipolygon" WKB definition. When creating polygons from points or lines, we desire to create a surface and this is why the corresponding geometry type is chosen wkbLinearRing for building polygons from both points or lines (rather than wkbPoint or wkbLine , respectively). The post-obit lawmaking block features an example for building a polygon shapefile delineating the hydraulic laboratory of the University of Stuttgart. The departure between the above example for creating a line shapefile are:
-
The project is
EPSG:4326. -
The point coordinates are generated within an
ogr.wkbLinearRingobject step-past-step rather than in a loop over dictionary entries. -
File, variable, and field names.
shp_dir = r "" + os . path . abspath ( '' ) + "/geodata/shapefiles/iws_va.shp" va_geo = create_shp ( shp_dir , layer_name = "basemap" , layer_type = "polygon" ) # create .prj file for the shapefile for GIS map applications with open ( shp_dir . split ( ".shp" )[ 0 ] + ".prj" , "w+" ) as prj : prj . write ( get_wkt ( 4326 )) # get basemap layer lyr = va_geo . GetLayer () # create polygon points pts = ogr . Geometry ( ogr . wkbLinearRing ) pts . AddPoint ( 9.103686 , 48.744251 ) pts . AddPoint ( 9.104689 , 48.744198 ) pts . AddPoint ( 9.104667 , 48.743960 ) pts . AddPoint ( 9.103557 , 48.744009 ) # create polygon geometry poly = ogr . Geometry ( ogr . wkbPolygon ) # build polygon geometry from points poly . AddGeometry ( pts ) # add together field to allocate building type field = ogr . FieldDefn ( "building" , ogr . OFTString ) lyr . CreateField ( field ) poly_feature_defn = lyr . GetLayerDefn () poly_feature = ogr . Characteristic ( poly_feature_defn ) poly_feature . SetGeometry ( poly ) poly_feature . SetField ( "edifice" , "Versuchsanstalt" ) lyr . CreateFeature ( poly_feature ) lyr = None va_geo = None
Build Shapefile from JSON¶
Loading geometry information from in-line defined variables is cumbersome in do, where geospatial information are often provided on public platforms (due east.g., land apply or cover). The following example uses a JSON file generated with map service data from the Baden-Württemberg State Institute for the Surroundings, Survey and Nature Conservation (LUBW), where polygon nodes are stored in WKB polygon geometry format ( "MultiPolygon (((node1_x node1_y, nodej_x, nodej_y, ... ...)))" ):
-
The JSON file (download hq100-dreisam.json and relieve it to a subdirectory called
/geodata/json/)is read with Pandas and the shapefile is created, as before, with thecreate_shpfunction. -
The projection is
EPSG:25832. -
2 fields are added in the course of
-
"tbg_name"(the original cord name of the polygons in the LUBW data), and -
"area"(a existent number field, in which the polygon area is calculated in m2 usingpolygon.GetArea()).
-
-
The polygon geometries are derived from the WKB*formatted definitions in the
"wkb_geom"field of the pandas dataframe objectdreisam_inundation.
# become data from json file dreisam_inundation = pd . read_json ( r "" + os . path . abspath ( '' ) + "/geodata/json/hq100-dreisam.json" ) # create shapefile shp_dir = r "" + os . path . abspath ( '' ) + "/geodata/shapefiles/dreisam_hq100.shp" dreisam_hq100 = create_shp ( shp_dir , layer_name = "basemap" , layer_type = "polygon" ) # create .prj file for the shapefile for GIS map applications with open ( shp_dir . split ( ".shp" )[ 0 ] + ".prj" , "due west+" ) every bit prj : prj . write ( get_wkt ( 25832 )) # get basemap layer lyr = dreisam_hq100 . GetLayer () # add cord field "tbg_name" lyr . CreateField ( ogr . FieldDefn ( "tbg_name" , ogr . OFTString )) # add string field "surface area" lyr . CreateField ( ogr . FieldDefn ( "area" , ogr . OFTReal )) for wkt , tbg in zilch ( dreisam_inundation [ "wkt_geom" ], dreisam_inundation [ "TBG_NAME" ]): # create Feature as child of the layer characteristic = ogr . Characteristic ( lyr . GetLayerDefn ()) # assign tbg_name feature . SetField ( "tbg_name" , tbg ) # use WKT format to add a polygon geometry to the Feature polygon = ogr . CreateGeometryFromWkt ( wkt ) # define default value of 0 to the area field feature . SetField ( "area" , polygon . GetArea ()) characteristic . SetGeometry ( polygon ) # append the new feature to the basement layer lyr . CreateFeature ( feature ) lyr = None dreisam_hq100 = None
Tip
Open the new dreisam_hq100.shp in QGIS and explore the attribute table.
Too GeoJSON data can be used to create an ogr.Geometry with ogr.createFromGeoJson(FILENAME) :
from osgeo import ogr geojson_data = """{"type":"Betoken","coordinates":[1013452.282805,6231540.674235]}""" betoken = ogr . CreateGeometryFromJson ( geojson_data ) print ( "Ten= %d , Y= %d (EPSG:3857)" % ( bespeak . GetX (), point . GetY ()))
X=1013452, Y=6231540 (EPSG:3857)
Calculate Geometric Attributes¶
The above code cake illustrates the usage of polygon.GetArea() to summate the polygon area in yard\(^2\). The ogr library provides many more functions to calculate geometric attributes of features and here is a summary:
-
Unify multiple polygons
wkt_... = ...
polygon_a = ogr.CreateGeometryFromWkt(wkt_1)
polygon_b = ogr.CreateGeometryFromWkt(wkt_2)
polygon_union = polygon_a.Marriage(polygon_b) -
Intersect 2 polygons
polygon_intersection = polygon_a.Intersection(polygon_b) -
Envelope (minimum and maximum extents) of a polygon
env = polygon.GetEnvelope()
impress("minX: %d, minY: %d, maxX: %d, maxY: %d" % (env[0],env[ii],env[1],env[3]) -
Convex hull (envelope surface) of multiple geometries (points, lines, polygons)
all_polygons = ogr.Geometry(ogr.wkbGeometryCollection)
for feature in POLYGON-SOURCE-LAYER: all_polygons.AddGeometry(characteristic.GetGeometryRef())
convexhull = all_polygons.ConvexHull()
Saveconvexhullto shapefile (use thecreate_shpfunction equally shown in the higher up examples or read more at pcjerick's Github pages)
Tip: To create a tight hull (e.thou., of a signal cloud), await forconcavehullfunctions. -
Length (of a line)
wkt = "LINESTRING (415128.5 5320979.5, 415128.6 5320974.5, 415129.75 5320974.7)"
line = ogr.CreateGeometryFromWkt(wkt)
print("Length = %d" % line.Length()) -
Area (of a polygon):
polygon.GetArea()(see above example) -
Example to calculate centroid coordinates of polygons.
Of import
The units of the geometric attribute (e.g., chiliad\(^two\), U.S. anxiety, or others) are calculated based on the definitions in the .prj file <prjp> (remember the definition of projections in WKT format in the Projections and Coordinate Systems section).
Export to Other Formats¶
The higher up examples bargain with .shp files merely, but other formats can exist useful (e.g., to create spider web applications or export to Google Globe). To this end, the following paragraphs illustrate the creation of GeoJSON and KML files. Several other conversions can be performed, not only between file formats but also between feature types. For instance, polygons can exist created from point clouds (amidst others with the ConvexHull method mentioned above). Interested students can larn more virtually conversions in Michael Diener's Python Geospatial Analysis Cookbook.
GeoJSON¶
GeoJSON files can be easily created as before, even without activating a commuter:
triangle = ogr . Geometry ( ogr . wkbLinearRing ) triangle . AddPoint ( - 11717151.498691 , 2356192.894805 ) triangle . AddPoint ( - 11717120.446149 , 2355586.175893 ) triangle . AddPoint ( - 11719392.059083 , 2354012.050842 ) polygon = ogr . Geometry ( ogr . wkbPolygon ) polygon . AddGeometry ( triangle ) with open up ( r "" + os . path . abspath ( '' ) + "/geodata/geojson/pitillal-triangle.geojson" , "w+" ) every bit gjson : gjson . write ( polygon . ExportToJson ())
For more than robust file treatment and defining a projection, activate the driver ogr.GetDriverByName("GeoJSON") . Thus, the creation and manipulation of GeoJSON files work similarly to the shapefile handlers shown above.
gjson_driver = ogr . GetDriverByName ( "GeoJSON" ) # make spatial reference sr = osr . SpatialReference () sr . ImportFromEPSG ( 3857 ) # create GeoJSON file gjson = gjson_driver . CreateDataSource ( "pitillal-full.geojson" ) gjson_lyr = gjson . CreateLayer ( "pitillal-full.geojson" , geom_type = ogr . wkbPolygon , srs = sr ) # get layer feature definitions feature_def = gjson_lyr . GetLayerDefn () # create new characteristic new_feature = ogr . Feature ( feature_def ) # assign the triangle from the above code cake new_feature . SetGeometry ( polygon ) # add together new feature to Layer gjson_lyr . CreateFeature ( new_feature ) # shut links to data sources gjson = None gjson_lyr = None
KML (Google Earth)¶
To display bespeak, line, or polygon features in Google Globe, features can be plugged into Google's KML (Keyhole Markup Language), similar to the creation of a GeoJSON file, and with the office geometry.ExportToKML :
triangle = ogr . Geometry ( ogr . wkbLinearRing ) triangle . AddPoint ( - 11717151.498691 , 2356192.894805 ) triangle . AddPoint ( - 11717120.446149 , 2355586.175893 ) triangle . AddPoint ( - 11719392.059083 , 2354012.050842 ) polygon = ogr . Geometry ( ogr . wkbPolygon ) polygon . AddGeometry ( triangle ) #geojson = poly.ExportToJson() with open ( r "" + os . path . abspath ( "" ) + "/geodata/pitillal-triangle.kml" , "w+" ) as gjson : gjson . write ( polygon . ExportToKML ())
Moreover, similar to GeoJSON files and shapefiles, KML files can be generated more robustly (with a defined project, for example). All you need to do is initiating the KML commuter ( kml_driver = ogr.GetDriverByName("KML") ) and ascertain a KML data source ( kml_file = kml_driver.CreateDataSource(FILENAME.KML) ).
Read Field Values of a Shapefile Gdal
Source: https://hydro-informatics.com/jupyter/geo-shp.html