-
Book Overview & Buying
-
Table Of Contents
-
Feedback & Rating

Geospatial Development By Example with Python
By :

Now that we have all we need installed, we will go through our first example. In this example, we will test the installation and then see a glimpse of OGR's capabilities.
To do this, we will open a vector file containing the boundaries of all the countries in the world and make a list of country names. The objective of this simple example is to present the logic behind OGR objects and functions and give an understanding of how geospatial files are represented. Here's how:
geopy
folder. Make sure that the data
folder is named data and that it's inside the geopy
folder.geopy
project opened, right-click on the project folder and select New | Directory. Name it Chapter1
.Chapter1
folder and select New | Python File. Name it world_borders.py
, and PyCharm will automatically open the file for editing.import ogr # Open the shapefile and get the first layer. datasource = ogr.Open("../data/world_borders_simple.shp") layer = datasource.GetLayerByIndex(0) print("Number of features: {}".format(layer.GetFeatureCount()))
world_borders
. An output console will open at the bottom of the screen, and if everything goes fine, you should see this output:C:\Python27\python.exe C:/geopy/world_borders.py Number of features: 246 Process finished with exit code 0
Congratulations! You successfully opened a Shapefile and counted the number of features inside it. Now, let's understand what this code does.
The first line imports the ogr
package. From this point on, all the functions are available as ogr.FunctionName()
. Note that ogr
doesn't follow the Python naming conventions for functions.
The line after the comment opens the OGR datasource (this opens the shapefile containing the data) and assigns the object to the datasource
variable. Note that the path, even on Windows, uses a forward slash (/
) and not a backslash.
The next line gets the first layer of the data source by its index (0
). Some data sources can have many layers, but this is not the case of a Shapefile, which has only one layer. So, when working with Shapefiles, we always know that the layer of interest is layer 0.
In the final line, the print statement prints the number of features returned by layer.GetFeatureCount()
. Here, we will use Python's string formatting, where the curly braces are replaced by the argument passed to format()
.
Now, perform the following steps:
# Inspect the fields available in the layer. feature_definition = layer.GetLayerDefn() for field_index in range(feature_definition.GetFieldCount()): field_definition = feature_definition.GetFieldDefn(field_index) print("\t{}\t{}\t{}".format(field_index, field_definition.GetTypeName(), field_definition.GetName()))
Number of features: 246 0 String FIPS 1 String ISO2 2 String ISO3 3 Integer UN 4 String NAME 5 Integer POP2005 6 Integer REGION 7 Integer SUBREGION Process finished with exit code 0
What happens in this piece of code is that feature_definition = layer.GetLayerDefn()
gets the object that contains the definition of the features. This object contains the definition for each field and the type of geometry.
In the for
loop, we will get each field definition and print its index, name, and type. Note that the object returned by layer.GetLayerDefn()
is not iterable, and we can't use for
directly with it. So first, we will get the number of fields and use it in the range()
function so that we can iterate through the indexes of the fields:
# Print a list of country names. layer.ResetReading() for feature in layer: print(feature.GetFieldAsString(4))
Number of features: 246 0 String FIPS 1 String ISO2 2 String ISO3 3 Integer UN 4 String NAME 5 Integer POP2005 6 Integer REGION 7 Integer SUBREGION Antigua and Barbuda Algeria Azerbaijan Albania Armenia Angola ... Saint Barthelemy Guernsey Jersey South Georgia South Sandwich Islands Taiwan Process finished with exit code 0
The layer is iterable, but first, we need to ensure that we are at the beginning of the layer list with layer.ResetReading()
(this is one of OGR's "gotcha" points).
The feature.GetFieldAsString(4)
method returns the value of field 4
as a Python string. There are two ways of knowing whether the country names are in field 4
:
Your complete code should look similar to the following:
import ogr # Open the shapefile and get the first layer. datasource = ogr.Open("../data/world_borders_simple.shp") layer = datasource.GetLayerByIndex(0) print("Number of features: {}".format(layer.GetFeatureCount())) # Inspect the fields available in the layer. feature_definition = layer.GetLayerDefn() for field_index in range(feature_definition.GetFieldCount()): field_definition = feature_definition.GetFieldDefn(field_index) print("\t{}\t{}\t{}".format(field_index, field_definition.GetTypeName(), field_definition.GetName())) # Print a list of country names. layer.ResetReading() for feature in layer: print(feature.GetFieldAsString(4))
Change the font size
Change margin width
Change background colour