Python 3 GIS | OGR | Fiona

Exercise 5: Working with vector data

You will need to download this data (http://drought.geo.msu.edu/teaching/geo429/data/mich_co_shp.zip)and unzip it to your folder to use it for this exercise.

Open and read the shapefile with OGR following the process we did in class. When using DumpReadable() method on a feature to print the feature into human readable format, Jupyter notebook doesn't show any output. This could be a problem with Jupyter notebook, or simply with the OGR python binding. So to work around this problem, you can try to run that piece of code within interactive python or python interpreter to get the list of attributes. Here is how to do it.

1. Open a Linux terminal on the server by clicking on the terminal icon
2. In the terminal, enter module load geo429 to set up the proper environment. Then enter python command to start the python interpreter. See the following example. Notice that you have to use Python 3.5.2 with Anaconda for ogr to work. Other versions might/might not work. Once you load the geo429 module, the python environment should be for 3.5.2.

   lluo@black:~$ module load geo429 lluo@black:~$ python Python 3.5.2 |Anaconda custom (64-bit)| (default, Jul 2 2016, 17:53:06) [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] on linux Type "help", "copyright", "credits" or "license" for more information. >>> 

3. Copy and paste your code into the Python interpreter to run the following code. It should now give you the output from DumpReadable method. The output gives you the attributes that are associated with this specific feature. In this shapefile, all geographic features are counties of Michigan, and each county has the same attributes. The actual values associated with these attributes will obviously vary. But this will give you an idea what attributes that you can extract and work with.

   >>>from osgeo import ogr >>>shpf = ogr.Open("data/mich_co.shp") >>>lyr = shpf.GetLayer(0) >>>feature = lyr.GetFeature(0) >>>feature.DumpReadable() OGRFeature(mich_co):0 FID (Integer64) = 0 w_M_201 (Real) = 5312.000000000000000 w_F_201 (Real) = 5176.000000000000000 b_M_201 (Real) = 19.000000000000000 b_F_201 (Real) = 27.000000000000000 nA_M_20 (Real) = 35.000000000000000 nA_F_20 (Real) = 41.000000000000000 a_M_201 (Real) = 6.000000000000000 a_F_201 (Real) = 19.000000000000000 pop (Real) = 10635.000000000000000 county (String) = Alcona marrigs (Real) = 51.000000000000000 divorcs (Real) = 33.000000000000000 POLYGON ((666667.467185152 480494.211763146,675142.281752073 480722.470884362,685792.431459327 481271.563931039,695290.517929201 481624.721626855,700062.989810847 481708.114326595,711456.115569082 482017.18594244,712100.355392683 480175.300987687,712974.328711306 479313.654057265,713359.156832663 477393.313839829,713885.474683868 476280.555491611,713672.894405104 474153.762459305,714188.458937549 471938.117692383,713694.791413845 470778.576325394,714144.89279161 468268.735545871,714774.799969475 467109.256103982,715988.963850573 466031.816244756,715674.904884499 463662.457350373,715212.114455894 461990.948710321,714620.334722815 458860.375639736,713424.491203851 456713.114505301,712899.162809242 454179.933320121,712873.64431283 452788.40893671,713210.606329985 448872.571019732,713693.291200278 446788.451035371,712947.017607005 443503.338960932,701836.417659922 443150.665716232,686638.68076447 442621.728127439,667811.141751631 441884.748812394,667480.119826987 453075.083321341,667076.870736058 464363.345038424,666841.937585521 475610.600646161,666667.467185152 480494.211763146))

   To exit from the Python interpreter, simply press Ctrl + d on the keyboard to exit. You can close the terminal once you have found out the above information.

4. Now you should have the information you need to work with the shapefile, i.e, attribute names.

Write code to find out the following using OGR:

• How many features are stored in this shapefile?
• Print out all the county name of each feature, the type of geometry for this county (i.e, point, line, polygon or multipolygon?) If you cannot do this with OGR, can you do it with Fiona?
• Write code to find the number of marriages and divorces in Ingham county.
• Write code to find the county with the highest population.

Working with Fiona

Using the Fiona module, open the mi_airports.shp file (http://drought.geo.msu.edu/teaching/geo429/data/mi_airports_shp.zip), and write code to find out the following.

• How many features are stored in this shapefile and what geometry type are they?
• How many airports are in Ingham county?
• Which airport has the highest elevation?

Improving matplotlib plotting

Use the functions that we developed in class, namely, getLineSegmentsFromGeometry, readShape, drawMultiPoints, drawMultiLines, drawMultiPolygons to work on this problem.

• Place all the functions (plotting and reading shapefiles with ogr) in a separate python file, called mymodule.py and save this mymodule.py in your current work directory. Congratulations! You just created your first module. You can also add DocString for your module to describe it. Your module will become very handy when we use these functions in the following questions. To use the function in your module, simply import them

  from mymodule import *

• Download this shapefile (http://drought.geo.msu.edu/teaching/geo429/data/nybb.zip) and unzip it. Import mymodule to read the shapefile. Make a plot of all the features in this shapefile using the functions in your module.
• Revise the plotting functions in your module, to allow users to specify the colors when drawing. Place your new plotting functions here, and replot the polygons in the nybb.shp shapefile.

Thanks in advance!