{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Working with vector data: OGR"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Vector data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sometimes, geospatial data is acquired and recorded for particular geometric objects such as polygons or lines. An example is a road layout, where each road is represented as a geometric object (a line, with points given in a geographical projection), with a number of added *features* associated with it, such as the road name, whether it is a toll road, or whether it is dual-carriageway, etc. This data is quite different to a raster, where the entire scene is tessellated into pixels, and each pixel holds a value (or an array of value in the case of multiband rasterfiles). \n",
"\n",
"If you are familiar with databases, vector files are effectively a database, where one of the fields is a geometry object (a line in our previous road example, or a polygon if you consider a cadastral system). We can thus select different records by writing queries on the features. Some of these queries might be spatial (e.g. check whether a point is inside a particular country polygon).\n",
"\n",
"The most common format for vector data is the **ESRI Shapfile**, which is a multifile format (i.e., several files are needed in order to access the data). We'll start by getting hold of a shapefile that contains the countries of the world as polygons, together with information on country name, capital name, population, etc. The file is available [here](http://aprsworld.net/gisdata/world/world.zip).\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will download the file with `wget` (or `curl` if you want to), and uncompress it using `unzip` in the shell:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--2014-11-11 16:00:23-- http://aprsworld.net/gisdata/world/world.zip\r\n",
"Resolving aprsworld.net... 72.251.203.219\r\n",
"Connecting to aprsworld.net|72.251.203.219|:80... connected.\r\n",
"HTTP request sent, awaiting response... 200 OK\r\n",
"Length: 3436277 (3.3M) [application/zip]\r\n",
"Saving to: `world.zip.3'\r\n",
"\r\n",
"100%[======================================>] 3,436,277 2.56M/s in 1.3s \r\n",
"\r\n",
"2014-11-11 16:00:25 (2.56 MB/s) - `world.zip.3' saved [3436277/3436277]\r\n",
"\r\n",
"Archive: world.zip\r\n",
" inflating: world.dbf \r\n",
" inflating: world.shp \r\n",
" inflating: world.shx \r\n"
]
}
],
"source": [
"# Downloads the data using wget\n",
"!wget http://aprsworld.net/gisdata/world/world.zip\n",
"# or if you want to use curl...\n",
"#! curl http://aprsworld.net/gisdata/world/world.zip -o world.zip\n",
"!unzip -o -x world.zip"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We need to import `ogr`, and then open the file. As with GDAL, we get a handler to the file, (`g` in this case). OGR files can have different layers, although Shapefiles only have one. We need to select the layer using `GetLayer(0)` (selecting the first layer)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from osgeo import ogr\n",
"\n",
"# have to make sure have access to gdal data files \n",
"import os\n",
"if 'GDAL_DATA' not in os.environ:\n",
" os.environ[\"GDAL_DATA\"] = '/opt/anaconda/share/gdal'\n",
"\n",
"g = ogr.Open( \"world.shp\" )\n",
"layer = g.GetLayer( 0 )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to see a field (the field `NAME`) we can loop over the features in the layer, and use the `GetField('NAME')` method. We'll only do ten features here:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"GUATEMALA\n",
"BOLIVIA\n",
"PARAGUAY\n",
"URUGUAY\n",
"SURINAME\n",
"FRENCH GUIANA\n",
"WESTERN SAHARA\n",
"GAMBIA\n",
"MOROCCO\n",
"MALI\n"
]
}
],
"source": [
"\n",
"n_feat = 0\n",
"for feat in layer:\n",
" \n",
" print feat.GetField('NAME')\n",
" \n",
" n_feat += 1\n",
" if n_feat == 10:\n",
" break"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you wanted to see the different layers, we could do this using:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Field 1: NAME\n",
"Field 2: CAPITAL\n",
"Field 3: APPROX\n",
"Field 4: AREA\n",
"Field 5: SOURCETHM\n"
]
}
],
"source": [
"layerDefinition = layer.GetLayerDefn()\n",
"\n",
"\n",
"for i in range(layerDefinition.GetFieldCount()):\n",
" print \"Field %d: %s\" % ( i+1, layerDefinition.GetFieldDefn(i).GetName() )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Each feature, in addition to the fields shown agove, will have a `Geometry` field. We get a handle to this using the `GetGeometryRef()` method. Geometries have many methods, such as `ExportToKML()` to export to KML (Google Maps/Earth format):"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"'-12.0443,14.669667 -11.87845,14.8252 -11.76455,15.039033 -11.63345,15.335633 -11.58515,15.634533 -11.52995,15.6713 -11.34705,15.6736 -11.26425,15.6529 -11.2125,15.549433 -11.1849,15.3092 -11.0503,15.097667 -10.9675,15.086167 -10.53275,15.365533 -9.70465,15.364367 -9.52175,15.373567 -9.4631,15.440233 -9.43205,15.619567 -9.37685,15.658667 -9.13185,15.656367 -9.1385,15.489433 -8.1165,15.487567 -6.83415,15.504167 -5.513,15.483867 -5.25265,16.2627 -5.2693,16.3033 -5.58505,16.439867 -5.63765,16.5635 -5.7623,17.728033 -5.8814,18.589933 -6.02265,19.525633 -6.1362,20.697567 -6.3107,21.8621 -6.4575,22.952833 -6.4686,23.543433 -6.52675,23.8627 -6.7123,24.997733 -5.05065,24.9834 -5.00065,24.982967 -4.51805,24.628633 -2.77295,23.4557 -2.17905,23.0869 -1.6281,22.7265 -0.5692,22.030967 0.6009,21.241133 1.0281,21.001933 1.10095,20.985233 1.21475,20.9807 1.21475,20.807767 1.24435,20.7486 1.40825,20.609033 1.5517,20.5514 1.6678,20.3724 1.9421,20.185067 2.0514,20.172933 2.2221,20.168367 2.27785,20.157 2.3348,20.1312 2.4031,20.028033 2.4486,20.009833 2.6307,19.971933 2.8208,19.887733 3.12355,19.778533 3.21345,19.6458 3.19755,19.4501 3.15885,19.3242 3.10535,19.174767 3.13495,19.065567 3.2738,18.875933 3.3091,18.853933 3.3774,18.852433 3.52535,18.9131 3.82355,18.996533 4.34825,19.149 4.34725,19.115667 4.3451,19.045533 4.29825,18.466367 4.23895,17.8664 4.1984,17.276833 4.1766,16.666467 4.1423,16.331633 3.97065,16.208933 3.8864,16.099767 3.85205,15.8315 3.8146,15.700467 3.6305,15.5809 3.4245,15.5081 3.07805,15.428033 2.75035,15.396833 2.3197,15.323 2.10435,15.2731 1.8297,15.216967 1.62995,15.213833 1.41145,15.120267 1.05565,15.0402 0.8809,15.046433 0.6905,15.0402 0.61245,14.965333 0.52195,14.890467 0.419,14.856167 0.238,14.858267 0.1196,14.8782 -0.2863,15.048533 -0.614,15.108833 -0.75445,15.1109 -0.82935,15.094267 -0.94485,15.018367 -1.285,14.7418 -1.4754,14.625333 -1.6439,14.5723 -1.94035,14.462067 -2.16195,14.312333 -2.24625,14.2094 -2.3118,14.0472 -2.42415,13.939067 -2.69255,13.824667 -2.77995,13.596967 -3.1014,13.524167 -3.3573,13.452433 -3.53205,13.248633 -3.6725,13.2341 -3.84105,13.232 -4.0564,13.2653 -4.1001,13.256967 -4.22805,13.0386 -4.24135,12.8486 -4.5398,12.277633 -4.72495,12.088433 -4.7833,12.04 -4.8209,12.017533 -4.8687,12.0 -5.03205,11.943633 -5.40485,11.807433 -5.4576,11.7805 -5.46435,11.758033 -5.39025,11.609867 -5.38015,11.431733 -5.43175,11.337433 -5.443,11.246133 -5.39135,11.051533 -5.40035,10.9767 -5.5171,10.836 -5.52835,10.806067 -5.5073,10.733367 -5.5285,10.681033 -5.5168,10.6259 -5.46805,10.575 -5.4553,10.5029 -5.46805,10.409567 -5.50265,10.376233 -5.5873,10.2946 -5.72715,10.208567 -5.8078,10.1817 -5.96645,10.1602 -6.0444,10.1602 -6.13045,10.187067 -6.1466,10.212167 -6.1466,10.260567 -6.0928,10.3466 -6.13855,10.4595 -6.17615,10.486367 -6.20305,10.486367 -6.238,10.477433 -6.2837,10.429033 -6.3469,10.4165 -6.40875,10.421867 -6.5405,10.4595 -6.5835,10.450533 -6.6319,10.423667 -6.6561,10.3878 -6.66415,10.307167 -6.7018,10.267733 -6.74755,10.256967 -6.8228,10.256967 -6.87255,10.239033 -6.9801,10.151233 -7.06885,10.1315 -7.20595,10.1315 -7.30005,10.1602 -7.33095,10.194267 -7.35785,10.310767 -7.41165,10.341233 -7.45735,10.353767 -7.7101,10.351967 -7.7531,10.4882 -7.8284,10.525833 -7.9427,10.529433 -8.19545,10.5366 -8.26805,10.545533 -8.33525,10.5814 -8.3339,10.619033 -8.28015,10.7176 -8.17795,10.921933 -8.20755,10.9524 -8.3124,10.982867 -8.385,10.9739 -8.42265,10.947033 -8.51945,10.807233 -8.627,10.771367 -8.67,10.7696 -8.68615,10.783933 -8.68615,10.8162 -8.6458,10.889667 -8.6458,10.904 -8.68885,10.9273 -8.6915,10.9775 -8.68615,11.038433 -8.666,11.102967 -8.6176,11.149567 -8.28955,11.3449 -8.4401,11.4435 -8.77625,11.5761 -8.86635,11.6263 -8.8986,11.674667 -8.88245,11.7123 -8.76415,11.8969 -8.748,11.998 -8.75615,12.0537 -8.9702,12.136333 -8.90365,12.335867 -9.1426,12.388267 -9.36645,12.279433 -9.7355,12.1323 -9.80655,12.0866 -9.96635,12.0 -10.0388,11.958267 -10.1407,11.925833 -10.31445,11.925867 -10.3944,11.966 -10.4546,11.986067 -10.4917,11.987633 -10.5241,11.975267 -10.5465,11.9559 -10.5583,11.948633 -10.56825,11.944533 -10.5777,11.943933 -10.58715,11.9458 -10.6018,11.9581 -10.6155,11.973867 -10.6342,12.0 -10.6706,12.034 -10.6969,12.058533 -10.72995,12.058533 -10.8209,12.019967 -10.8457,12.019967 -10.8953,12.053 -10.92835,12.149433 -10.91185,12.3257 -10.9201,12.3615 -10.97385,12.358733 -11.12265,12.237567 -11.17225,12.033733 -11.4244,12.141167 -11.56495,12.254067 -11.60215,12.303667 -11.3913,12.429 -11.4946,12.489567 -11.4326,12.610767 -11.5401,12.6824 -11.5649,12.7237 -11.57315,12.756767 -11.4905,12.988133 -11.5194,13.096933 -11.7385,13.336567 -11.8625,13.325567 -11.9989,13.344833 -12.08155,13.416433 -12.08985,13.4798 -12.0774,13.866767 -12.0443,14.669667'"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"the_geometry = feat.GetGeometryRef()\n",
"the_geometry.ExportToKML()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Many of the methods that don't start with `__` are interesting. Let's see what these are. typically, the interesting methods start with an upper case letter, so we'll only show those:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"AddGeometry\n",
"AddGeometryDirectly\n",
"AddPoint\n",
"AddPoint_2D\n",
"Area\n",
"AssignSpatialReference\n",
"Boundary\n",
"Buffer\n",
"Centroid\n",
"Clone\n",
"CloseRings\n",
"Contains\n",
"ConvexHull\n",
"Crosses\n",
"Destroy\n",
"Difference\n",
"Disjoint\n",
"Distance\n",
"Empty\n",
"Equal\n",
"Equals\n",
"ExportToGML\n",
"ExportToJson\n",
"ExportToKML\n",
"ExportToWkb\n",
"ExportToWkt\n",
"FlattenTo2D\n",
"GetArea\n",
"GetBoundary\n",
"GetCoordinateDimension\n",
"GetDimension\n",
"GetEnvelope\n",
"GetEnvelope3D\n",
"GetGeometryCount\n",
"GetGeometryName\n",
"GetGeometryRef\n",
"GetGeometryType\n",
"GetPoint\n",
"GetPointCount\n",
"GetPoint_2D\n",
"GetPoints\n",
"GetSpatialReference\n",
"GetX\n",
"GetY\n",
"GetZ\n",
"Intersect\n",
"Intersection\n",
"Intersects\n",
"IsEmpty\n",
"IsRing\n",
"IsSimple\n",
"IsValid\n",
"Length\n",
"Overlaps\n",
"PointOnSurface\n",
"Segmentize\n",
"SetCoordinateDimension\n",
"SetPoint\n",
"SetPoint_2D\n",
"Simplify\n",
"SimplifyPreserveTopology\n",
"SymDifference\n",
"SymmetricDifference\n",
"Touches\n",
"Transform\n",
"TransformTo\n",
"Union\n",
"UnionCascaded\n",
"Within\n",
"WkbSize\n"
]
}
],
"source": [
"for m in dir ( the_geometry ):\n",
" if m[0].isupper():\n",
" print m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You'll notice that many of these mechanisms e.g. `Overlaps` or `Touches` are effectively geoprocessing operations (they operate on geometries and return `True` if one geometry overlaps or touches, respectively, the other). Other operations, such as `Buffer` return a buffered version of the same geometry. This allows you to actually do fairly complicated geoprocessing operations with OGR. However, if you want to do geoprocessing in earnest, you should really be using [Shapely](http://toblerity.org/shapely/manual.html). \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A particularly useful webpage for this section is [available in the OGR cookbook](http://pcjericks.github.io/py-gdalogr-cookbook/). Have a look through that if you want more in depth information."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Selecting attributes and/or data extents"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"OGR provides an easy way to select attributes on a given layer. This is done using a SQL-like syntax (you can read more on [OGR's SQL subset here](http://www.gdal.org/ogr/ogr_sql.html). The main point is that the *attribute filter* is applied to a complete layer. For example, let's say that we want to select only countries with a population (field APPROX) larger than 90 000 000 inhabitants:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PAKISTAN has 123490000 inhabitants\n",
"JAPAN has 124710000 inhabitants\n",
"RUSSIAN FEDERATION has 150500000 inhabitants\n",
"INDIA has 873850000 inhabitants\n",
"BANGLADESH has 120850000 inhabitants\n",
"BRAZIL has 159630000 inhabitants\n",
"NIGERIA has 91700000 inhabitants\n",
"CHINA has 1179030000 inhabitants\n",
"INDONESIA has 186180000 inhabitants\n",
"JOHNSTON ATOLL has 256420000 inhabitants\n",
"KINGMAN REEF - PALMYRA ATOLL has 256420000 inhabitants\n",
"UNITED STATES has 256420000 inhabitants\n"
]
}
],
"source": [
"g = ogr.Open ( \"world.shp\" )\n",
"lyr = g.GetLayer( 0 )\n",
"lyr.SetAttributeFilter ( \"APPROX > 90000000\" )\n",
"for feat in lyr:\n",
" print feat.GetFieldAsString ( \"NAME\") + \" has %d inhabitants\" % \\\n",
" feat.GetFieldAsInteger(\"APPROX\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So we get a list of popoulous countries (note that Johnston Atoll and Palmyra are part of the US, and report the sample popuation as the US!)\n",
"\n",
"An additional way to filter the data is by geographical extent. Let's say we wanted a list of all the countries in (broadly speaking) Europe, *i.e.* a geographical extent in longitude from 14W to 37E, and in latitude from 72N to 38N. We can use `SetSpatialFilterRect` to do this:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ALGERIA ---- ALGIERS\n",
"BELGIUM ---- BRUSSELS\n",
"LUXEMBOURG ---- LUXEMBOURG\n",
"SAN MARINO ---- SAN MARINO\n",
"AUSTRIA ---- VIENNA\n",
"CZECH REPUBLIC ---- PRAGUE\n",
"SLOVENIA ---- LJUBLJANA\n",
"HUNGARY ---- BUDAPEST\n",
"SLOVAKIA ---- BRATISLAVA\n",
"YUGOSLAVIA ---- BELGRADE [BEOGRADE]\n",
"BOSNIA AND HERZEGOVINA ---- SARAJEVO\n",
"ALBANIA ---- TIRANE\n",
"MACEDONIA, THE FORMER YUGOSLAV REPUBLIC ---- SKOPJE\n",
"LITHUANIA ---- VILNIUS\n",
"LATVIA ---- RIGA\n",
"BULGARIA ---- SOFIA\n",
"BELARUS ---- MINSK\n",
"MOLDOVA, REPUBLIC OF ---- KISHINEV\n",
"IRELAND ---- DUBLIN\n",
"ICELAND ---- REYKJAVIK\n",
"SPAIN ---- MADRID\n",
"SWEDEN ---- STOCKHOLM\n",
"FINLAND ---- HELSINKI\n",
"TURKEY ---- ANKARA\n",
"RUSSIAN FEDERATION ---- MOSCOW\n",
"GREECE ---- ATHENS\n",
"PORTUGAL ---- LISBON\n",
"POLAND ---- WARSAW\n",
"NORWAY ---- OSLO\n",
"GERMANY ---- BERLIN\n",
"ESTONIA ---- TALLINN\n",
"TUNISIA ---- TUNIS\n",
"CROATIA ---- ZAGREB\n",
"ROMANIA ---- BUCURESTI\n",
"UKRAINE ---- KIEV\n",
"NETHERLANDS ---- AMSTERDAM\n",
"JERSEY ---- SAINT HELIER\n",
"GUERNSEY ---- SAINT PETER PORT\n",
"FAROE ISLANDS ---- TORSHAVN\n",
"DENMARK ---- COPENHAGEN\n",
"MONACO ---- MONACO\n",
"ANDORRA ---- ANDORRA LA VELLA\n",
"LIECHTENSTEIN ---- VADUZ\n",
"SWITZERLAND ---- BERN\n",
"ISLE OF MAN ---- DOUGLAS\n",
"UNITED KINGDOM ---- LONDON\n",
"FRANCE ---- PARIS\n",
"VATICAN CITY (HOLY SEE) ---- VATICAN CITY\n",
"ITALY ---- ROME\n"
]
}
],
"source": [
"g = ogr.Open ( \"world.shp\" )\n",
"lyr = g.GetLayer( 0 )\n",
"lyr.SetSpatialFilterRect ( -14, 37, 38, 72)\n",
"for feat in lyr:\n",
" print feat.GetFieldAsString ( \"NAME\") + \" ---- \" + feat.GetFieldAsString ( \"CAPITAL\") "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Saving a vector file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Saving a vector file using OGR requires a number of steps:\n",
"\n",
"1. Definition of the format\n",
"2. Definition of the layer projection and geometry type (e.g. lines, polygons...)\n",
"3. Definition of the data type of the different fields\n",
"4. Creation of a feature, population of the different fields, and setting a geometry\n",
"5. Addition of the feature to the layer\n",
"6. Destruction of the feature\n",
"\n",
"This appears quite involved, but let's see how this works. Note that when you generate a new vector file, OGR will fail if the file already exists. You might want to use `os.remove()` to get rid of the file if it exists. \n",
"\n",
"Let's see how this is done with an example which is a snippet that creates a GeoJSON file with the location of the different national parks. GeoJSON is a nice geographic format, and [github allows you to display it easily as a map](https://github.com/blog/1528-there-s-a-map-for-that)."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# National park information, separated by TABs\n",
"import os\n",
"from osgeo import ogr,osr\n",
"\n",
"\n",
"parks = \"\"\"Dartmoor national park\\t-3.904\\t50.58\n",
"New forest national park\\t-1.595\\t50.86\n",
"Exmoor national park\\t-3.651\\t51.14\n",
"Pembrokeshire coast national park\\t-4.694\\t51.64\n",
"Brecon beacons national park\\t-3.432\\t51.88\n",
"Pembrokeshire coast national park\\t-4.79\\t51.99\n",
"Norfolk and suffolk broads\\t1.569\\t52.62\n",
"Snowdonia national park\\t-3.898\\t52.9\n",
"Peak district national park\\t-1.802\\t53.3\n",
"Yorkshire dales national park\\t-2.157\\t54.23\n",
"North yorkshire moors national park\\t-0.8855\\t54.37\n",
"Lake district national park\\t-3.084\\t54.47\n",
"Galloway forest park\\t-4.171\\t54.87\n",
"Northumberland national park\\t-2.228\\t55.28\n",
"Loch lomond and the trossachs national park\\t-4.593\\t56.24\n",
"Tay forest park\\t-4.025\\t56.59\n",
"Cairngorms national park\\t-3.545\\t57.08\"\"\"\n",
"\n",
"# See if the file exists from a previous run of this snippet\n",
"if os.path.exists ( \"parks.json\"):\n",
" # It does exist, so remove it\n",
" os.remove ( \"parks.json\" )\n",
"\n",
"# We need the output projection to bet set to Lat/Long\n",
"latlong = osr.SpatialReference()\n",
"latlong.ImportFromEPSG( 4326 )\n",
"\n",
"# Invoke the GeoJSON driver\n",
"drv = ogr.GetDriverByName( 'GeoJSON' ) \n",
"# This is the output filename\n",
"dst_ds = drv.CreateDataSource( 'parks.json' )\n",
"# This is a single layer dataset. The layer needs to be of points\n",
"# and needs to have the WGS84 projection, which we defined above\n",
"dst_layer = dst_ds.CreateLayer('', srs =latlong , \\\n",
" geom_type=ogr.wkbPoint ) \n",
"\n",
"# We just need a field with the Park's name, and its type is a String\n",
"field_defn=ogr.FieldDefn( 'name', ogr.OFTString )\n",
"dst_layer.CreateField( field_defn )\n",
"\n",
"\n",
"# Algorithm is as follows:\n",
"# 1. Loop over lines\n",
"# 2. Split line into park name, longitude, latitude\n",
"# 3. Create WKT of the point\n",
"# 4. Set the attribute name to name of park\n",
"# 5. Clean up\n",
"\n",
"for park_id, line in enumerate( parks.split( \"\\n\" ) ):\n",
" # Get the relevant information\n",
" park_name, lon, lat = line.split(\"\\t\")\n",
" # Create a geogrpahical representation of the current park\n",
" wkt = \"POINT ( %f %f )\" % ( float(lon), float(lat) )\n",
" # Create a feature, using the attributes/fields that are\n",
" # required for this layer\n",
" feat = ogr.Feature(feature_def=dst_layer.GetLayerDefn())\n",
" # Feed the WKT into a geometry\n",
" p = ogr.CreateGeometryFromWkt( wkt )\n",
" # Feed the geometry into a WKT\n",
" feat.SetGeometryDirectly( p )\n",
" # Set the name field to its value\n",
" feat.SetField ( \"name\", park_name )\n",
" # Attach the feature to the layer\n",
" dst_layer.CreateFeature( feat )\n",
" # Clean up\n",
" feat.Destroy()\n",
"\n",
"# Close file \n",
"dst_ds = None\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see the result of this on [github](https://gist.github.com/jgomezdans/6811102). \n",
"\n",
"Additionally, note that if we had defined a coordinate transformation as in the raster session, we could apply this transformation to an OGR geometry entity (in the snippet above, `p` would be such), and it would be reprojected.\n",
"\n",
"**Exercise**\n",
"Modify the above snippet to output a GeoJSON file for the Peak District National Park, whose UTM30N ([EPSG code: 32630](http://spatialreference.org/ref/epsg/32630/)) co-ordinates are $577659, 5911841$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Rasterising"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A very frequent problem one finds is how to mask out an area in a raster file that is defined as polygon in a shapefile. For example, if you have a raster of the worlds population density, and you want to extract all the pixels that belong to one particular country, how do you go about that? One way around this is to *rasterise* the polygon(s), which translates into \"burning\" pixels that fall within the polygon with a number, resulting in a mask.\n",
"\n",
"The way to do this is to use GDAL's `RasterizeLayer` method. The method takes a handle to a GDAL dataset (one that you create yourself, with the right projection and geotransform, as you've seen above), and a OGR layer. The syntax for `RasterizeLayer` is\n",
"\n",
" err = gdal.RasterizeLayer ( raster_ds, [raster_band_no], ogr_layer, burn_values=[burn_val] )\n",
" \n",
"where `raster_ds` is the GDAL raster datasource (note that it needs to be georreferenced, *i.e.* it requires projection and geotransform), `raster_band_no` is the band of the GDAL dataset where we want to burn pixels, `ogr_layer` is the vector layer object, and `burn_val` is the value that we want to burn.\n",
"\n",
"Let's use `gdal.RasterizeLayer` in conjunction with all that we have covered above. Say we want to create a mask that only selects the UK or Ireland in `world.shp`, and we want to apply this mask to the MODIS landcover product that we used in the GDAL session (h17v03 tile ), file `lc_h17v03.tif`. We find that in this case, `world.shp` is in longitude latitude, and the MODIS data is in the MODIS projection, so we will reproject the vector data to match the MODIS data (so the latter is not interpolated and artifacts introduced). To make this efficient and avoid saving to disk, we shall use *in-memory vector and rasters*, and we will output a numpy array as our mask. Note then the steps:\n",
"\n",
"1. Crate the projection conversion object (as for GDAL before)\n",
"2. Create an in memory **raster** dataset to store the mask, using `lc_h17v03.tif` as a reference for geotransforms, array size and projection.\n",
"3. Create an in memory **vector** dataset to hold the features that will be reprojected\n",
"4. Open `world.shp` and apply an `AttributeFilter` to select a country\n",
"5. Select a geometry from `world.shp`, project it and store it in the destination in memory vector layer\n",
"6. Once this is done, use `gdal.RasterizeLayer` with both in-memory raster and vector datasets\n",
"7. Read the in memory raster into an array\n",
"\n",
"This is a particularly good exercise that will stress all that we have learned so far."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQkAAAEACAYAAACgZ4OsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGOpJREFUeJzt3U+MJOV9xvHvzw4otmFEUJJuYFYBKevIezJCmTk4lq0o\nIXCZwReDpVjIRp5IxH+UHALkYiu5OJYcxT4YaWRsg5Wsg2SZGRRwAMuRLEWZCRGYddbr3ZW8ETMw\njUGOxj5lSX45dNVMbVFdXV1dXfVW1fORVvRUd0+/3cz79Puv3jJ3R0Rkkrc1XQARCZtCQkRyKSRE\nJJdCQkRyKSREJJdCQkRy1RoSZnaHmZ0zswtm9kCdry0i5Vhd6yTM7O3AT4A/APaBfwc+4u4/rqUA\nIlJKnS2JFeCiu19y98vAt4D1Gl9fREqoMyRuAl5O/LwXHRORgNUZElr/LdJCv1Lja+0DJxI/n2Dc\nmjhiZgoSkYa4u2UdrzMkngdOmtnNwCvA3cBH0g9aWVmpsUjV2dvbY3l5uelizKyt5Yb2ln2R5R4M\nBoxGo5mft7u7O/G+2kLC3d80s08C/wy8HXhEMxsi1SoTENPU2ZLA3Z8Gnq7zNUVkPlpxWZGlpaWm\ni1BKW8sN7S1728qtkKhI2/7Hx9pabmhv2dtWboWEiORSSIhILoWEiORSSIi02GAwWPhrKCREWmwR\n6yLSFBIikkshISK5FBIigapjvKEIhYRIg/KCoI7xhiIUEiINygqCwWAQTCsCaj7BS0Qmi4MhlBZE\nTCEh0pBkKJTdB6IO6m6I1CzuToQaCmkKCZEaJccakkER/ze08QhQSIjULt2CSAbDaDQKroWhkBBp\nUKiDlUkKCZEapcMg+XNo3YyYZjdEapLXpQh5IFMtCZGaJAcnJ90XIoWESINC7WIkKSREapTVYgi5\nFQEKCZHaJUMh9IAAhYSITKHZDZGatG3AMqaQEMmRrNhlKnQb9ouYRiEhUlC6wrdhIVQVFBIiE0yr\n+PMEQ1taEaCBS+m4EM+qhHa1PBQS0lnp07KlHIWEdFLe+EGZ51etTd0NjUlIp6jFUD2FhHTGpIAI\n6Vs7pLIUpZCQTqgiINTFyKYxCWm1IrMXRSv/IitxWwMCFBLSYiFU/ja8/rwUEtJZIVTOEMowr7lC\nwswumdlLZvaCme1Gx643s2fN7LyZPWNm1yUe/5CZXTCzc2Z2+7yFF8kz6xhDlWMSIe56Xda8LQkH\nPujut7r7SnTsQeBZd3838L3oZ8zsFHA3cAq4A/iKmaklI6UUqdBNDVp2JRxiVVRSS/28Bjwa3X4U\nuCu6vQ6cdvfL7n4JuAisIFLCtIo4a0WNv/nnreBdCwiopiXxnJk9b2afiI4N3D3+pEZAHNE3AnuJ\n5+4BN835+tJjkyp1UxW1iwEB86+TeJ+7v2pmvwE8a2bnkne6u5uZ5zz/Lfft7R3nyNLSEktLS3MW\nUbps3ovuxs9J/rcPDg8POTw8LPTYuULC3V+N/vszM/sO4+7DyMyG7n5gZjcAr0UP3wdOJJ6+HB27\nwvLy8jxFkh6ZdF3Nos9Nb3E/76nf6cAJWfoL+JVXXpn42NLdDTN7p5ldG91+F3A7cAbYBu6NHnYv\n8ER0exu4x8yuNrNbgJPAbtnXl35L7xhV9a5RZX9XGwJiVvO0JAbAd8ws/j1/7+7PmNnzwONmdh9w\nCfgwgLufNbPHgbPAm8D97p7XFRGZaJ5v7EV2KboWEAAWUj01M19Z0YSHTDZLMOQ9tuqgyAutvH0y\nQ2l57O7u4u7pmUpAJ3hJC6QrdJErceeFQFUBkRzoTFf2Sa+RdXyW8jQRKFrMJEGb9+zORW1jn54J\nqWtT3CZmXxQSErRFrJqsoqKlAyI+G7WLU6gKCQnWpG7FpOAoss5hUWMRXaaQkGDN2orIW++wqIrc\nREDU/ZoKCWmVabMH054XwkzCvOp+DwoJCd60EOhCxQ+ZQkI6oWhLIZR1CWVpClRkDnkLp4qsrQhd\nU2XXYippjSoqSVtnIpoMN7UkpHOKLMVuc4uibgoJ6bSssYo+rG2okrob0kmTQkEtiNmpJSHBq6pi\nJ5dPS3EKCQlamRZAVhAkux1ta000XV51NyRYi/jGb1sroumAALUkJHBlK0nZ5dvzvm4XKSQkWPMG\nxDythra1OBZJISFBqqKSFj3FPFShlFchIcFZ1PkVah2Uo5CQ4FQdEG0Lh9BmYDS7IZ00afNcmZ1a\nEtIZWRfsKfKNHNK3dojUkpDWy1pwVXQRVmjncYQYWGpJSKvlBUT6dlpoAREqtSSktYrOgmStmwgx\nIEJsRYBCQlqm6CXz0lfXSj8npIAINRxi6m5Iq+TtZVlk8VSye9Jk5ZxlYLVpCglpnSKDkTD9eqBN\ntCaSZQupNZNHISGdNGm8oulv7rYEQ5LGJKQz8q4+ngyNkCpqG3bMUkhIZ6TDIH1fyEK+Hoi6G9Ip\nbb04T4hliqklIb0RB8Sk7kadaydCDoU0tSREIgqIbAoJ6Zysyp63CKvolGoV2hYQoO6G9EB6liOp\nSBej7N6YRa4Y1obZjaktCTP7mpmNzOxM4tj1ZvasmZ03s2fM7LrEfQ+Z2QUzO2dmtyeO32ZmZ6L7\nvlT9WxEZm+Vs0Kq7GPEiqXj8I29VZdET0ZpWpLvxdeCO1LEHgWfd/d3A96KfMbNTwN3Aqeg5XzEz\ni57zMHCfu58ETppZ+neK1GpS5a1qyXQyMNIrLCd1iUIMi6kh4e4/AH6eOrwGPBrdfhS4K7q9Dpx2\n98vufgm4CKya2Q3Ate6+Gz3uscRzRCqVPrErfcm/or9jUdf9KNK9CSkwyg5cDtw9/rRHQPxubgT2\nEo/bA27KOL4fHRepVNaZn1n3J6Wv7hVK5YQwuiFzz264uwNeQVlESpv0zZvXckivmwihQoao7OzG\nyMyG7n4QdSVei47vAycSj1tm3ILYj24nj+9n/eK9veMGx9LSEktLSyWLKH02y3hCuvUQYmui6tmP\nw8NDDg8PCz22bEtiG7g3un0v8ETi+D1mdrWZ3QKcBHbd/QA4NLPVaCDzo4nnXGF5efnonwJCipi2\nLmKSSWMVIQXEoiwtLV1R1/IUmQI9Dfwr8Dtm9rKZfQz4PPCHZnYe+P3oZ9z9LPA4cBZ4Grg/6o4A\n3A98FbgAXHT375Z6dyIZps1GTLrKeJG1DCFociDTjutw88zMV1ZWmi6GtNS0Cp/XpcjaCi9Eiwqy\n3d1d3N2y7tOKS+mUqipRiHtPQDMrNHXuhnRG0YozrZUhV1JIiLRQnaGmkBAh7EHLpikkJChbW1sL\n+b2TFloVPQGrzxQSEpyqg6JMZVdAHFNISFDW19cZDocLWxdQpFuhgLiSQkIaEbcWkv/d2to6qqDb\n29tsb29X9nptuVpWiBQSUqtkV2JnZ+eo1QAwHA7fEgw7OztzX4lL4TAfLaaS2mxtbbGxscFgMGA4\nHLK2tnZ033A4zHxO8jFJky7EA/OFQpu6GnVdHkAhIQuVbDnELYW44m9vb3NwcMDGxgYHBweZQVFm\n5WPZytOmgID6WkgKCVmorIqf7FLE92c9bnV1FWhf5a2LWhLSSltbW6yvrx/dnibdnYh/nuePvw+t\niDopJGRhJo0zpCW7ILF5Km2ol/JbhDreq2Y3pFLr6+sMBgN2dnYKPb7Kac4yQtpwNlQKCalc/M02\naWaiDrNW/Kqv0tWldRkKCVmI1dXVwq2ERYXJrDMibbXosiskpHKzdDcWbdaFVnkX7Mn7HVmthzYH\nT5JCQiq1tbU1UzO7ji5J0aBIX38j/Tu61IWYhWY3pBLJlkMorYhY2Yo9aWv9vgWFQkLmVubU7ipb\nEHnXyZhlS7uqphOb6GYscipU3Q2Zy2AwYGNjY+L9WWFQVUBMa/7HFX+eac42zVQsKpy0pb5UIrmL\nc/KPdXNzM3NRVRVBMet5HXVV9KYHLMu8T22pLwtXdJYAqg2IWfRlJWbV2+4rJKRy8R/n1tbWFa2I\nKrsZsaKX86vr273pVkRSVafPKyRkIZIBkXXS1rznZsz6+NAuAly3eVpRCgmpTDzLMenErj5U0pDD\nqGxQaHZDKjMcDqfuH1FEvL9lVSd/zVtxmzwPpGplPgeFhByZtN6hyDqIaX98ZSt8FUFRdiBPZ4iO\nqbvRY3Hlj0/vHg6HU1dLJjeVSdrc3Jz6eln7RmQ9RhZr1m6HQqLH4s1oiyyjjv+wkkEyaV/KEJVp\nRcz7el1phWgxVU+lpyenOTg4AIrvNjVNVosirxUx6/RplTtml/1doYdE8n3lLabSmESPra6uFq58\nkwYly0oHwrRuRhu7IWVOO69T0RBTd6On4vMt4nGCRVbCSd2SeV6zyPhGFaresSo0RYJC3Y2eSi4w\nauO3dFIcQlXstF3VKsX078vanyIUo9FI3Q15q3g2IsRvt1kNh8OjMZN530/VZ31O2sAmFEXep7ob\nPZMcsAzpj7WM5JXANjY2OhF4IVJLomfW19c5ODi4ooI1uav1POJv6dXV1dYERGjBXKQ8U0PCzL5m\nZiMzO5M49jkz2zOzF6J/dybue8jMLpjZOTO7PXH8NjM7E933pRLvRyqyvr7emkrVJaEFRFFFWhJf\nB+5IHXPgb9391ujf0wBmdgq4GzgVPecrZhYPhjwM3OfuJ4GTZpb+ndKA+Ft4dXX16NqbbaFl0/WY\nGhLu/gPg5xl3ZY2ErgOn3f2yu18CLgKrZnYDcK2770aPewy4q1yRpSohDajFA4+zakuLKOQwm1a2\necYkPmVmPzSzR8zsuujYjcBe4jF7wE0Zx/ej49Kgra0ttra22NnZOVpq3dT4RJmFWm2aum1LmGUp\nO7vxMPBX0e2/Br4I3FdFgfb2jrNkaWmJpaWlKn6tTJCsnG2qdG0U0vkcFy5c4PDwkGuuuWbqY0uF\nhLu/Ft82s68CT0Y/7gMnEg9dZtyC2I9uJ4/vZ/3u5eXlrMPSkDadxBW6UAIC4OTJk1f8fP78+YmP\nLdXdiMYYYh8C4pmPbeAeM7vazG4BTgK77n4AHJrZajSQ+VHgiTKvLeUl94UoskdEGwIitAsBTRJS\nQCQVWTg2tSVhZqeBDwC/bmYvA58FPmhm72U8y/FT4E8A3P2smT0OnAXeBO7343Xf9wPfAN4BPOXu\n3y3zpmQ+Rc/+rCog2hA0faZzN+QKXThPI62K8zXqEGpLIvbkk0/q3A0ptntU7ODgoPS0pHSLQqJn\n4oo/baqzTV2E0FsRbafuRo9kNXnr2E9i0Ra1UrSq08ZD72pAfndDZ4H2SPIPPevqWm0Lirg1tIjL\n91VVsdsQENMoJHoqfSXwvIAIeYaizosAT1rGPunKZPH0YvpY24JDYxI9NcvmKqEGxPb2duUVLu+k\nsbzjWc/LenzbAgLUkui9nZ0d1tbWWF1dPZoiTc5qxNfY6OL0aZY4NGcNiknaGAppComeSw76xaeM\nSzu7BYuikOipaZe+29nZYXV1tRXLnhcxcAnHQdH3wFBI9FRepYrP62hDQCxa3wMCNHApkbgizHpl\nr6QurtDse0CAQkIio9ForoCAcGdBZD4KCQEmXy28z9SKGFNICHA81Vm2y9ClroY22L2SBi7lyDwD\nlU12Naqc2ci7onhfg0MhIa22qIDQmaXH1N0QoJ3TnVXv7F3lNUC7RCHRc/GW+pO09RKA85gUFH0N\nEIVEj02b8lxbW2M0GgUZFHHZZPEUEj1VZE1EshKura0FGRZ10sCl9EK85Ho4HB59GyfP49jZ2eHg\n4ICNjY2gxynUiqiPWhI9ErceNjY2ji4UnDwFfDAYsLa2xnA4DP608L5+qzdBe1z2WMgthWmaOqW9\nq+GkLfXlLRQQUpTGJHqkC7tLNTl42tVWxDRqSfRIqNOZRWnasxkKiR5qa1DMcgUyqY5CoqPiqc6s\nq4e3vcvRhL52NUBjEp0UT3XGax7i/0LYG8MkLxCUdbGgpsre54AAhURnxJvGpFdSxrdDDodYPN4Q\nd4fioGiqe9T3cIgpJDqkjdOacQAkWw3JwUmth2ieQqJj4ovsxOIKN+/+lfNKtgayZijiMFAohEch\n0RHJ/SnTFwaGZrsboU9dKiDyKSQ6KDmjkRcOdV0IOB5XCC0oFA7FKCQ6ZpZuRV2ti9ACQuEwG4WE\nlJKecUhPV66trbG5uclwOAwmIBQO5WgxVYfs7OwUah1UMaUYh8KkhVmj0Yj19fVgTsZSQJSX25Iw\nsxPAY8BvAg5suvuXzex64B+B3wIuAR929/+OnvMQ8HHgf4FPu/sz0fHbgG8Avwo85e6fWcQb6qus\nlZVpWdON84jHGpLLpUMJhbS4NaOwmF3ufhJmNgSG7v6imV0D/AdwF/Ax4HV3/4KZPQD8mrs/aGan\ngH8Afhe4CXgOOOnubma7wCfdfdfMngK+7O7fTb2e9pOYw9bWFhsbG0c7TMWSLYeql2Qnp1xD6Vbk\nUUhky9tPIrcl4e4HwEF0+5dm9mPGlX8N+ED0sEeBfwEeBNaB0+5+GbhkZheBVTP7L+Bad9+NnvMY\n47C5IiRkPslp0DgYkt+gVQdE+jVCF2pAhH7l8sIDl2Z2M3ArsAMM3D3+yxgB8Tu8Efi3xNP2GIfK\n5eh2bD86LgsyLRzKdj2mLYpqq/R7qbPShhwQUDAkoq7Gt4HPuPsvzI5bJVFXIpw98ORI3OWYVLHL\nLuNuazjM8o0dP3aR3/KhtyBiU0PCzK5iHBDfdPcnosMjMxu6+4GZ3QC8Fh3fB04knr7MuAWxH91O\nHt/Per29veMGx9LSEktLSwXfiiTFA5lZAVEmHEIdkJxVenfwoo9dZDma8Prrr/PGG28Ueuy0gUtj\nPObwhrv/WeL4F6Jjf2NmDwLXpQYuVzgeuPztqLWxA3wa2AX+CQ1cLkT87Tfpvlm7F/HAZFtbD/NY\n1MWDmw6ILKUHLoH3AX8MvGRmL0THHgI+DzxuZvcRTYECuPtZM3scOAu8Cdzvxyl0P+Mp0HcwngLV\noOUC5FXmzc3NzH0akuIWQ3JKtY8BAdW/77Z+jtpSvyfib67Nzc2jadL4WKjnVoSuSGugyQHRWczT\nkpCOGI1GRxvTwJXdkq6MN9RtWrehLQExjZZl98ik08mlvKKfY1sDAtSSEJnbtKBoc0CAWhIiMoVC\nQmSB2t6KAIWEyEwGg8HRv+SxLtOYhPTetAVoec8r8ri2U0tCeq9MQPSJQkJkAk0TjykkRCZQS2JM\nISESSQ5IKiCOKSREUhQQV1JIiEguhYRIS9U1sKqQEGmpujYD0mIqEdo5DhEHRJFT0ucJE4WE9F7b\nAiJrG73kseRO6VW0NBQS0muLDIisClpmt+4ivyMvLOalkJDeWkRAFNmBu8hr57UCsnbEWuTYhEJC\npAKzVtJZWwqTXq+OrpJCQnqlidZD3vOygmKW31fH7IamQKU3QhygbMNJZAoJEcmlkBCZQxtaAvNS\nSEjnpbebq/p3d51CQiQyGo2O/sU/i2Y3pONmucJW1n15FwnuS4ioJSGdVTYgqnxOF6glIZ0xy5Ln\nsvoYFGpJVOTw8LDpIpTS1nLDcdmLDEyGNM7Qts9cIVGRtv2Pj7W13PDWsicHHbOEEBDQvs9cISGd\nMy0sZDYak5DWUyAslrl702U4YmbhFEakZ9zdso4HFRIiEh6NSYhILoWEiOQKJiTM7A4zO2dmF8zs\ngabLk2Zml8zsJTN7wcx2o2PXm9mzZnbezJ4xs+sSj38oei/nzOz2Gsv5NTMbmdmZxLGZy2lmt5nZ\nmei+LzVY9s+Z2V70ub9gZneGVnYzO2Fm3zez/zSzH5nZp6Pjrfjcp3L3xv8BbwcuAjcDVwEvAu9p\nulypMv4UuD517AvAX0S3HwA+H90+Fb2Hq6L3dBF4W03lfD9wK3CmZDnjcapdYCW6/RRwR0Nl/yzw\n5xmPDabswBB4b3T7GuAnwHva8rlP+xdKS2IFuOjul9z9MvAtYL3hMmVJj/6uAY9Gtx8F7opurwOn\n3f2yu19i/EewUkcB3f0HwM/nKOeqmd0AXOvuu9HjHks8Z2EmlB3e+rlDQGV39wN3fzG6/Uvgx8BN\ntORznyaUkLgJeDnx8150LCQOPGdmz5vZJ6JjA3ePJ+lHQLw2+EbG7yHW9PuZtZzp4/s0W/5PmdkP\nzeyRRJM9yLKb2c2MW0M7tP9zB8IJiTbMw77P3W8F7gT+1Mzen7zTx+3DvPcRxHssUM7QPAzcArwX\neBX4YrPFmczMrgG+DXzG3X+RvK+Fn/uRUEJiHziR+PkEVyZq49z91ei/PwO+w7j7MDKzIUDUVHwt\nenj6/SxHx5oySzn3ouPLqeONlN/dX/MI8FWOu21Bld3MrmIcEN909yeiw6393JNCCYnngZNmdrOZ\nXQ3cDWw3XKYjZvZOM7s2uv0u4HbgDOMy3hs97F4g/uPYBu4xs6vN7BbgJOMBqabMVE53PwAOzWzV\nzAz4aOI5tYoqV+xDjD93CKjs0es8Apx1979L3NXaz/0KTY+cJkaI72Q8KnwReKjp8qTKdgvj0egX\ngR/F5QOuB54DzgPPANclnvOX0Xs5B/xRjWU9DbwC/A/jcZ6PlSkncBvjCnkR+HJDZf8448G7l4Af\nMq4wg9DKDvwe8H/R38cL0b872vK5T/unZdkikiuU7oaIBEohISK5FBIikkshISK5FBIikkshISK5\nFBIikkshISK5/h9DvU8RB+NrbQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from osgeo import ogr,osr\n",
"import gdal\n",
"\n",
"reference_filename = \"lc_h17v03.tif\"\n",
"target_vector_file = \"world.shp\"\n",
"attribute_filter = \"NAME = 'IRELAND'\" \n",
"burn_value = 1\n",
"\n",
"# First, open the file that we'll be taking as a reference\n",
"# We will need to gleam the size in pixels, as well as projection\n",
"# and geotransform.\n",
"\n",
"g = gdal.Open( reference_filename )\n",
"\n",
"# We now create an in-memory raster, with the appropriate dimensions\n",
"drv = gdal.GetDriverByName('MEM')\n",
"target_ds = drv.Create('', g.RasterXSize, g.RasterXSize, 1, gdal.GDT_Byte)\n",
"target_ds.SetGeoTransform( g.GetGeoTransform() )\n",
"\n",
"# We set up a transform object as we saw in the previous notebook.\n",
"# This goes from WGS84 to the projection in the reference datasets\n",
"\n",
"wgs84 = osr.SpatialReference( ) # Define a SpatialReference object\n",
"wgs84.ImportFromEPSG( 4326 ) # And set it to WGS84 using the EPSG code\n",
"\n",
"# Now for the target projection, Ordnance Survey's British National Grid\n",
"to_proj = osr.SpatialReference() # define the SpatialReference object\n",
"# In this case, we get the projection from a Proj4 string\n",
"\n",
"# or, if using the proj4 representation\n",
"to_proj.ImportFromWkt( g.GetProjectionRef() )\n",
"target_ds.SetProjection ( to_proj.ExportToWkt() )\n",
"# Now, we define a coordinate transformtion object, *from* wgs84 *to* OSNG\n",
"tx = osr.CoordinateTransformation( wgs84, to_proj )\n",
"\n",
"# We define an output in-memory OGR dataset\n",
"# You could also do select a driver for an eg \"ESRI Shapefile\" here\n",
"# and give it a sexier name than out!\n",
"\n",
"drv = ogr.GetDriverByName( 'Memory' ) \n",
"dst_ds = drv.CreateDataSource( 'out' )\n",
"# This is a single layer dataset. The layer needs to be of polygons\n",
"# and needs to have the target files' projection\n",
"dst_layer = dst_ds.CreateLayer('', srs = to_proj, geom_type=ogr.wkbPolygon ) \n",
"\n",
"# Open the original shapefile, get the first layer, and filter by attribute\n",
"vector_ds = ogr.Open( target_vector_file )\n",
"lyr = vector_ds.GetLayer ( 0 )\n",
"lyr.SetAttributeFilter( attribute_filter )\n",
"\n",
"\n",
"# Get a field definition from the original vector file. \n",
"# We don't need much more detail here\n",
"feature = lyr.GetFeature(0)\n",
"field = feature.GetFieldDefnRef( 0 )\n",
"# Apply the field definition from the original to the output\n",
"dst_layer.CreateField( field )\n",
"feature_defn = dst_layer.GetLayerDefn()\n",
"# Reset the original layer so we can read all features\n",
"lyr.ResetReading()\n",
"for feat in lyr:\n",
" # For each feature, get the geometry\n",
" geom = feat.GetGeometryRef()\n",
" # transform it to the reference projection\n",
" geom.Transform ( tx )\n",
" # Create an output feature\n",
" out_geom = ogr.Feature ( feature_defn )\n",
" # Set the geometry to be the reprojected/transformed geometry\n",
" out_geom.SetGeometry ( geom )\n",
" # Add the feature with its geometry to the output yaer\n",
" dst_layer.CreateFeature(out_geom )\n",
" # Clear things up\n",
" out_geom.Destroy\n",
" geom.Destroy\n",
"# Done adding geometries\n",
"# Reset the output layer to the 0th geometry\n",
"dst_layer.ResetReading()\n",
"\n",
"# Now, we rastertize the output vector in-memory file\n",
"# into the in-memory output raster file\n",
"\n",
"err = gdal.RasterizeLayer(target_ds, [1], dst_layer,\n",
" burn_values=[burn_value])\n",
"if err != 0:\n",
" print(\"error:\", err)\n",
"\n",
"# Read the data from the raster, this is your mask\n",
"data = target_ds.ReadAsArray()\n",
"\n",
"\n",
"\n",
"# Plotting to see whether this makes sense.\n",
"\n",
"ndata = g.ReadAsArray()\n",
"plt.imshow ( ndata, interpolation='nearest', cmap=plt.cm.gray, vmin=0, vmax=1, alpha=0.3 )\n",
"plt.hold ( True )\n",
"\n",
"plt.imshow ( data, interpolation='nearest', cmap=plt.cm.gray, alpha=0.7 )\n",
"plt.grid ( False )\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Using matplotlib to plot geometries"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using matplotlib to plot geometries from OGR can be quite tedious. Here's an example of plotting a map of Angola from the `world.shp`. In the same vein of recommending Shapely and Fiona above for serious geoprocessing of vector data, you are encouraged to use [descartes](https://bitbucket.org/sgillies/descartes/) for plotting vector data!"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAM0AAAEACAYAAAAQvfK3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl4VEW6/z9v9oSEJSBESJAtiASQIEscZNMLGMVf8DoI\nCAqOdxAXFgURGBVwRJFcBtQr13EwMoMC4oIOw7AJsonsYZPdsIctSBaSTjrprt8fCbkBOyHd6e7T\n3anP85wn3Wep9+2c/nadU6fqW6KUQqPRVB4/oxPQaLwNLRqNxk60aDQaO9Gi0WjsRItGo7ETLRqN\nxk5cIhoRSRaRQyKyV0S+EZFaroij0RiBq2qa1UCcUupu4CgwyUVxNBq34xLRKKXWKKWsJW+3AdGu\niKPRGIE77mn+APzbDXE0GrcQ4OiBIrIGiLKxabJSalnJPn8CzEqphY7G0Wg8DXFV3zMRGQ78EXhA\nKZVfzj6645vGY1FKia31rmo9exB4BUgqTzBlEqv0MmXKFLv2d2TxhRi+8BmMjlERrrqn+QAIB9aI\nSKqIzHVRHI3G7Th8T1MRSqlYV5Sr0XgCXtUjoGfPnjqGB5Rf3WO4rCGgUsFFlJHxNZryEBGUOxsC\nNBpfRotGo7ETLRqNxk60aDQaO9Gi0WjsRItGo7ETLRqNxk60aDQaO9Gi0WjsRItGo7ETLRqNxk60\naDQaO9Gi0WjsRItGo7ETl4pGRMaJiFVEIl0ZR6NxJy4TjYjEAL2BU66KodEYgStrmr8AE1xYvkZj\nCK5yo0kCziql9rmifI3GSFxhFvgnir2b+5Tdvbxypk6dWvq6Z8+ebhkXrtHczPr161m/fn2l9nW6\nR4CItAHWAnklq6KBc0BnpdSlm/bVHgEaj6QijwCXG2uIyAngHqXUrza2adFoPBKjjTW0KjQ+hbZw\n0mhsYHRN4xYyMjLo27cvw4cPNzoVjY/jM6L56KOPsFqtrF69mrVr1xqdjsaHcYmXs7vIz89n+fLl\nrF+/noULFzJ79myysrJ4/PHH+eijjxgwYIDRKWp8EK++p+nTpw9Xr16lffv2/P73v6devXoAHD58\nmAkTJtCtWzeGDBlCUFAQffr0QaTcx0UazQ0Y2uRcEVUVTWxsLG+//TbNmjX7zbZ169YxYcIEevXq\nxRdffMFtt91WlVQ11YyKROOVl2dWa/EcuGFhYZjN5hu2ZWVlMWfOHFJTU1m+fDkPPfSQESlqfBiv\nawhYvHgxUVFR3HXXXRw9erT0kkspxcqVKxk0aBCNGzfm4MGDWjAal+BVl2cFBQWEhITwxhtv8Msv\nv9CjRw86dOhARkYGb731FlevXuWTTz6hS5cuLsxaUx3wmec03333HQBvvvkmn3/+OTVr1mTz5s0M\nHTqUnj17kpqaqgWjcTledU+zefNmAFq0aEHnzp2Jjo5m0KBBLFiwgKFDhxqcnaa64FWXZ/n5+YSG\nhvLUU08xevRolFIsXryYzz77jLS0NMLCwlyYraY64TOXZ4WFhQClzcciwuDBg8nMzDQyLU01w6tE\nExERwbhx4zh16v9sB6xWK2azmZCQEAMz01QnvEo0ABMnTmTt2rWcO3cOALPZTHBwMH5+XvdRNF6K\n133T6tWrx4svvsi8efOA/7vP0WjchdeJ5tNPP+W9995j3bp1AOTl5ekGAI1bcVmTs4iMAp4HLMBy\npdSrzij3yy+/ZODAgSQkJACwd+9eGjVq5IyiNZpK4RLRiEgv4P8B7ZRShSLitN6S2dnZ7Nq1i4sX\nL/Lpp59y5MgRvvrqK2cVr9HcEpc8pxGRJcBHSql1t9jP7l7Ohw8fJjU1FZPJRK1atejbty/h4eGV\nPt5qtZKRkUF6ejomk4mwsLDSxWq1YjKZyMvLIy8vr9zXAQEBREZG3rDUqVOH4OBg8vPzb1gKCgp+\ns+7mxWQyceeddzJo0CA9fMFDcPvQABFJBb4DHgTygfFKqZ029nOZR8Dx48fZvn07hw4d4tChQ5w+\nfZrz589z6dIlatSoQYMGDQgODqagoACTyYTJZMLPz4+QkJDfLMHBwTcsFouFnJwccnJyyM7OJisr\ni6ysrNKWvOtLUFBQ6XLz+8DAwBteb9y4ET8/PxYsWECnTp1c8j/RVB6XDA24hVlgAFBHKZUgIp2A\nJcBvB73gHLNApRT79+/n8OHDBAcHk5OTw8svv0xcXBwtWrSgffv2JCYmUr9+ferWrUtwcLDdMVzN\nsGHD6N69O0eOHNGiMQBDzQIBRGQFMEMptaHk/XGgi1Lqyk37OVTTHDx4kMjISKKioli2bBkTJ04k\nKyuL1q1bYzabCQgIIDEx0evcOp9++mnCwsIYN26cHqptMEZcnj0LNFRKTRGRlsD3SqnGNvazWzQr\nVqxgwIABiAiBgYH4+/szZcoUEhISvP5+4PTp0+zatYu3336bwsJC/P39jU6p2mLEyM0UIEVE9gNm\n4ClnFRwVFUVwcDCTJk2iffv21KhRg6CgIGcVbyiNGzcmODiYevXqacF4MF7Vy/k6W7duJSkpiUmT\nJtGtWzcXZGYcBQUF3H///WRmZur+dAbiM72cr5OQkMDXX3/NzJkzyc/PNzodp7J27VqaN2+u+9J5\nMF57Zu677z5+97vfsWjRIqNTcSqFhYWEhob6zCWnL+KVl2fXOX78OJ07d+azzz6jQYMGTszMOHJy\ncujXrx/Z2dm6tjEQn7s8u06LFi148cUXmT59OleuXLn1AV5AeHg4d9xxB8uWLTM6FU05eLVoAF57\n7bXSLijXjTe8GRFhxIgRPPvssxw5csTodDQ28HrRBAUFkZKSwpYtW0hJSWHFihVGp1Rlrt+veUtt\nU1hYyPr161m8eDGbNm3i6tWrRqfkUrzKjaYiWrVqxcqVK+nZsydRUVHEx8cbnVKVMJlMXLt2zeg0\nAEhPT2fnzp3k5uZy2223ER4eXtqXzmq1Mn36dHbu3Enz5s25fPkyqampfPnll/z+9783OnWX4DOi\nAWjTpg2ff/45Q4cOZf78+URF2eoa5x0MGDCAiRMn0qtXL3r06GFYHtOnT2fGjBmlD5IzMzMxmUyY\nzeZSo5P27duTkpJCjRo1ABgyZIjPNMzYwqtbz8pj8uTJpKWl8corrzi9bHcyc+ZMunbtyujRow3L\n4eGHH6ZHjx488MADldr/8uXLDB48mMuXLxMQ4L2/yT7belYeTZo0oaCgwOg0qkzjxo356aefDM3h\n+lCIynLu3Dlatmzp1YK5FT75yc6cOeMTU2skJiYyePBglixZwuOPP25IDs2aNePMmTOV3t9qtXL8\n+HFeeOEFgoKCCAgIQERKB/GZTCbCw8Px9/fHYrFgsVjw8/Ojc+fOJCYmEh0d7cJP4xx8sqZp2rQp\nR48eNTqNKlOrVi3effddXn31VR577DEuXLjg9hzatWvHyZMnK71/fHw8EyZMICwsDKVU6cjUoKAg\noqKiuPPOO4mMjCQiIoK6desSFRVFZGQk33zzDW3btiUuLo60tDTXfSAn4JP3NCaTiZiYGFJSUnzC\ndCM/P5/k5GQCAwP58ssv3Rr7X//6F++++y5z5sxxeSyLxcL7779PWFgYH3/8scvjVUS1u6cJDQ1l\nyJAhPmO4ERISQrdu3Th69Cgmk8ltcc1mM59//jmxsbFuiefv78/w4cNZtGgReXl5bonpCD4pGoBX\nXnmFf/7znz7TvaZr166EhYW55Rf/Oo888gjp6elunZGhTp061KhRg4sXL7otpr34rGiio6N58skn\nWbhwodGpOIXAwEDi4uLc9mU6c+YMO3bsYObMmdSqVcstMQGOHDmCn5+fRzcIuEQ0ItJZRLaLSKqI\n7Cgx13A7I0eOZPXq1Rh53+ZMkpKS+Mc//sHZs2ddHmvz5s106NDB7U3HX331FS+88AKBgYFujWsP\nrqppZgKvK6XigTdK3rud1q1bEx4ezqFDh4wI73QaNmxI7969mT9/vstjnTx50u2NKJmZmaxbt44R\nI0a4Na69uEo054HrdXpt4JyL4lSIiNCtWzcOHjxoRHiX0KNHD/7nf/7H4SZ1k8nEoUOHbln7nj17\nlnr16jkUw1Hee+89hg4d6vHP2FwlmonALBE5DSQDk1wU55bEx8dz+PBho8I7nS5duvDAAw8wfvx4\nu1vS9u/fz2233Ub37t259957SU9Pt7mf2WwmNTXVrTVNRkYG69atY8aMGW6L6SiuMgscDYxWSi0V\nkQEUu9P0tlWOM8wCK6JXr17MmDEDpZTXWzxd5/nnn+ell16iTZs2vPrqq3To0IH69esTGRlJcHAw\nVquVrKwsjh07xqpVq9i5cyeHDh3i8uXLPPfccwwYMID58+cTFxdHp06diIuLKzXxuHLlChs2bKBB\ngwZu7Sleu3ZtgoODuXTpkl02w87CE8wCs5VSNUteC5CplPpNE4wrbWmvo5SiYcOGfPDBBzRt2tSl\nsdyJ1Wpl69atLFu2jPT0dH799VcyMzMxm82ICDVr1qR+/fokJCRw11130bx5c8LDw6ldu3bpzf2F\nCxc4cuQIJ0+eLO1fFh4eTvPmzenQoYPbf2T++te/snv3blavXk1kZKRbY9+MEWaBu4GXlFIbROQB\nit02f9OC5g7RAEyZMoVdu3Yxbdo0l8fyBLy1VlVK8f7773Pw4EHWrVtHzZo1DcvFCNF0BD4EggET\n8LxSKtXGfm4RTU5ODi1btuSdd96hbdu2Lo+ncRylFG+++SbNmzdn9uzZhuXh9m40SqmdSqkuSqn2\nSql7bQnGnURERJCcnMzMmTPt6uaucT8iwsMPP8zWrVuNTqVcfLZHQFlycnJ45plnOHToEKmphupX\nUwmio6M9uqdztRBNSEgI//Vf/wUUt9JoPJvs7Gzq1q1rdBrlUi1EExgYyIcffsj777/PrFmzjE5H\ncwvOnz9P48a/mWTCY6gWorlOUlISBw4c8Jm+aL7K+fPnadKkidFplEu1Ec2cOXNo3749Tz/9tFc2\nx1YnTp06RbNmNifO8wh8cuTmdSwWC9u2bSM2NpamTZvy+eefe3SXcw1s2bKFt956i+3bt3PHHXcY\nlke1G7l5nWnTptG1a1fuvPNOevXqpQXjBSxdupRZs2YZKphb4bOiMZlMvPPOOwDk5uYa6h2mqTxZ\nWVkeLRjwUQsnoNT98Y477mDMmDGG92XSVI6AgACPseMtD5+taWrWrMn//u//EhwcTNeuXY1OR1NJ\nevbs6fF9BH1WNADPPPMMZrOZU6dOGZ2KppJcuXLFkKEB9uCzl2dQ3AISGxtLenq6Rzdhaoq5evUq\nX3/9NQcOHDA6lQrx6ZoGoEGDBoY4U2rs5+LFiwQFBbnV/cYRfF40Q4cOJSUlxS0OLpqq0apVK7p0\n6cJrr71mdCoV4vOi6du3L9OmTWPs2LEUFRUZnY7mFjz11FN8++23RqdRIQ6LRkQGiMjPImIRkQ43\nbZskIsdE5LCI9Kl6mlVjxIgRFBQUeHR3c00xNWvWJDs726P7B1alptkPPApsLLtSRFoDA4HWwIPA\nXBExtEb7+9//DkDz5s2NTENTCWrXrk39+vWZOnUq+fn5RqdjE4e/zEqpw0opW+ZbScAipVShUuok\ncBzo7GgcZ5CYmIjJZPIZX2dfxt/fn1mzZvHDDz/Qrl07j7TfckUN0BAoe9d9FjB0vouoqCifmkXA\n14mKimL27NkMHDiQHj16cOLECaNTuoEKRSMia0Rkv43lETvjGH6BGhcX5zP2tNWF/v378+ijj3rc\n3KkVPtxUStk0+LsF54CYMu+jqcCW1tVmgdf529/+xmOPPeaSsjWuY8iQIfTp0wez2UxQUJDL4rjV\nLFBEfgDGK6V2lbxvDSyk+D6mEfA90MLWwBl3WTgBJCcn87e//Y3x48dzzz33uCWmxjk8/vjjfP31\n17Rv395tMV0ynkZEHhWRM0ACsFxEVgAopQ4CS4CDwAqKPc8MvzwbP348EyZM4M9//jMLFy706CZN\nzY0kJCR41DxDPj1y0xZpaWkkJibSpk0bJk2apIc+ewHnz5/nqaee4pdffqFOnTpuiVltR27aolmz\nZuzevZuNGzfq3s9ewu23307btm1Zs2aN0akA1VA0UGzplJub6/b5VzSOk5mZ6THD1aulaCwWCyEh\nISxdulT3R/MCCgoKOHHiBHfeeafRqQDVVDShoaFs27aNn3/+mSeeeIJz5wyZqE1TSf71r3/RuXNn\nj3HdrJaiAYiNjWX16tU88sgjfPzxx0ano6kAq9VKVJSt+cOModqK5joLFy70mGtljW02bNhAUlKS\n0WmUUq1FIyL069ePFStW6Es0D0UpxZ49e3jwwQeNTqWUai0agJSUFMaNG8eYMWN0o4AHcuXKFWrU\nqEFERIRDxx8/fpxXX32Vjh07Eh8fz4QJEzhz5gyXL19m+vTpTJs2jby8PLvKrPaiAXj66ae5cOEC\nubm5RqeiuYnc3NzfuNPk5uby66+/lr63Wq039PBQSrFjxw6eeOIJunTpwsWLF3n22Wd58cUXOX/+\nPG3btqV58+akpqby7bffsmDBArty8mk3msqydetW2rRp4/GGDtWRRo0acfHiRXJzc9m0aRPTpk1j\n7969+Pv7Ex8fj1KKLVu2EBAQQGxsLK1ateLo0aNkZmbSv39/vvnmmxtE1759e4YPH05AQADh4eGM\nGzfO7py0aICmTZty6tQprFYrfn668vUkAgICaNasGQ8//DBpaWmMGjWKWbNmISJs2LCBgIAAZs6c\nidVqJS0tjdOnT9O7d286dOhQ7rksO7FXWlqa3T3rtWiAFi1aEBYWxokTJ/SQaA/kpZdeYtu2bUyZ\nMuWGGZ97975x5Err1q1p3bq1XWX7+fnZ3f9Qi4Zik7orV67QsGFDo1PR2CA+Pp74+HiXlW+vF4G+\nFgGOHj1KdHQ0oaGhRqeicTOdO3dm+fLldh2jRQO0a9eO8+fPa+ONakhBQYHdww20aACz2Wx0ChoD\nuHbtGhs2bKBfv352HVcl0dxkGHhPmfW9RWSniOwr+durKnFczYYNG2jdurXHdAjUuJ7MzEzGjh3L\noEGD7O5GVdWGgOuGgX/lRseZy0A/pdQFEYkDVlFssOGRFBYWEhISYnQaGjeydetW0tPTmTVrlt3H\nVqmmKc8wUCm1Ryl13ar/IBAqIoFVieVK7rrrLo4fP250Gho30rNnT5o0aUJKSordx7rjnuYxYJdS\nqtANsRxi/vz5NG7c2Og0NG5CKcVnn33G0aNHHXoud8vLMxFZA9gazDBZKbXsFsfGATMAR/zT3MLJ\nkydZsGCBQ9W0xnv56KOPSEtLo2nTpnYfe0vROGgYiIhEA98ATyqlyvUVdZdZoC2ys7NJSEhgyJAh\nHjOUVuN6RAQRueHqwq1mgSVJ3GwYWBvYAExRSpU72YgRFk7Xyc3NZciQIWRnZ5OcnGxIDhrjeOih\nh9i+fXu506+7zMKpPMNA4EWgOTBFRFJLFo+yfpk3bx5XrlzhzTffNDoVjQH07NmTBx98kO+++85u\n48hqZxZYVFTEnDlzSE5OZtq0aXTq1Mmt8TWegVKKzZs389577zFq1ChefvnlG7ZXVNNUuw6bhw4d\n4q233mLOnDm0bdvW6HQ0BiEidOrUiYYNG5KdnW3XsdWuG02NGjWoUaOGFkw159KlS4wYMYKYmBgm\nTpxo17HVpqb56aefWLNmDfv379fPZDTMnDmThx9+mHfffdfu8TTV4p5m8eLFjB49mri4OGJiYhgx\nYoQeBlCNMZlM9O7dm6tXr5bbfapa39Ps3r2bMWPGkJycTJs2bYxOR+MB/PLLL7Ro0cLh/oY+e09j\nMpkYOXIkffr0YdSoUVowmlLq16/P2bNnOXDggEPH+6xo5s+fz+7du1m8eDGJiYlGp6PxIOrXr8/4\n8ePp0aMHM2bMsPs5jc9enlmtVu644w63TQKkcS6FhYVs2rSp1NPsZm8zWwQHB9+wBAUFUVBQQF5e\nHnl5eeTm5mIymTCZTJjNZuLj45k0aRJPPvkkjRpVfgJynxSNUoqtW7d6lGm2xj5+/PFH/vznP/Mf\n//EfpX3FKnKOsVgsFBQUlIoiPz+f/Px8QkNDiYiIICIigvDwcGrWrElYWBh169YlOjqaCRMm2CUY\n8EHRHDlyhD/84Q+6T5mXIyLcd999fPPNN0an8ht8SjSnT59m7NixZGZm8umnn+Lv7290ShoHKWny\nNToNm/hMQ8CuXbu4++678ff3Z9KkSVowXo4ni8ZnaprMzEyCg4N57bXXCAz02JHVmkriyaLxmZrG\n39+fzMxMCgs9dlS1xg6Kioo89sfPJ0Rz6tQp+vfvzzvvvENYWJjR6WicQH5+vseeS58QzYkTJ4iN\njaV79+5Gp6JxEvn5+dSoUcPoNGzisGhuMgrsYGN7YxG5JiL2TwBiJxcvXiQyMtLVYTRupKioiIAA\nz7zlrkpNc90ocGM52/8C2Ocs7SC1atWyeyCRxrMJCQmxe1o/d+GwlJVShwGbT2hFpD+QBrhlPr7s\n7OzfTDGn8W5CQkIwmUxGp2ETp9/TiEg4MAGY6uyybbFjxw727t3rsb9KGscICgryWNFUWNM4aBQ4\nFZitlMqTSgyJq4rv2aJFi3j55ZeJiYmx2/ld49lYrVa3PqB2q+9ZiefZOKXU7pL3G4GYks21ASvw\nulJqro1jqzRyc/DgwbRo0YL+/fs7XIbGM1m7di1btmzh22/Ltc1zKe4YuVlauFKqtN1XRKYAObYE\n4wxMJpO+l/FRLBaLx3aFqkqTc3lGgW7j0qVLN0xcqvEdLBaLxzY5V6X1bCmw9Bb7THO0/FuRkZHB\nvn377J7NV+Md+GRNYzShoaFYLBbtKuOjuLshwB68VjTBwcGIiMc2S2qqRlFREUFBQUanYROvFc3q\n1atp2bKlbgjwUcxms8deRXitaBYvXkyfPn2MTkPjIgoKCrRonI2IeOx4C03V8dUOm4biyf9UjW/j\nlaJZtWoVy5cvJy4uzuhUNC7Ck4c7e91Ptdls5pFHHmH27NkOzcyr0VQVr6tpzpw5Q7169UhISDA6\nFY0L0U3OTmTJkiVaMNWAwsJCgoODjU7DJl4nmu3bt9OxY0ej09C4GF3TOIm0tDQ2bNigGwA0huJV\nopk7dy5JSUlER0cbnYrGxfj5+WGxWIxOwyZeJZpjx47RqlUro9PQuAEtGidw6dIldu7cqSeZrSb4\n+flhtVqNTsMmXiGajRs30qZNG/r160fLli2NTkfjBvz8/CgqKjI6DZu4xCxQRNqJyE8ickBE9omI\nw22HSimGDRvG5MmTGTFihKPFaLwMT748q0qPgOtmgX8tu1JEAoAFwFCl1H4RqQPY7UqulOKTTz5h\nw4YNpRP8aKoPSqlyZz0zGleYBfYB9iml9pfsd9WR8s+ePcsf//hHAgMDSU5O9th/oMY1mEwmj/V/\ncEXfs1hAichK4DZgsVLK7nn8br/9dpo1a8aIESN0LVMNKSgoICQkxOg0bOIKs8BA4D6gI2AC1orI\nLqXUOls7l2cWGBAQwMKFC+nXrx9t27alQYMGt/goGl/C39/frQ0BRpsFDgQSlVLDS96/BuQrpf7b\nxrG3NAucOnUqGzZs4L//+zeHa3wUq9XKuHHjGDZsGM8884whOVRkFuisJueyha8C2opIaEmjQA/g\nZ0cLnjRpEjt27CAnJ6eqOWq8hBUrVrBlyxa6detmdCo2cbpZoFIqk+JpNnYAqcAupZTDRoInT55E\nRCgoKHC0CI2X8dBDDzFgwACeffZZj5wOssqXZ1UKXonLs86dO9OzZ08GDhzopqw0noDFYmHcuHG0\na9eODz/80O3x3XF55jKio6O5du2a0Wlo3Iy/vz9vvfUWy5YtY9my8tqcjMEjRfP6669Tq1YtVqxY\nweTJk1m1apXRKWkMIDw8nDFjxjB16lSP8gvwSNFs27aN6OhoVq5cicVi8dhZfjWup1u3buTk5DB3\nrksmnnAIjxRNo0aN+PXXX4FiAcXGxhqckcYo/P39SU5O5s0332TRokVGpwN4aEPApk2b6N69O82b\nN8disTB27FjdK6Cac/z4ccaOHctrr73GCy+84PJ4FTUEeKRolFIMHjyYL774grFjxzJ06NBKl2mx\nWFi5ciVpaWmEhYURHx9PfHy87rvmA5w7d47Ro0czcOBA3n77bZfOKuB1ooFi4cTExPD666/bZaQx\ne/ZsDh8+zH/+53+SmZnJt99+S3h4OEOGDGHfvn3k5OQwfPhwGjVq5KyPoXEjV69eZfLkydSuXZvF\nixcTGRnpkjhe2eS8d+9e8vPzufvuuyt9zKZNm9i4cSPff/89kydPZubMmRw6dIhXXnmFRYsWYbFY\naNmyJS+99JILM9e4kjp16vDBBx8QFRVFhw4d2Lt3r9tz8NiaJi8vj9jYWKZPn07btm1vWVZGRgZP\nPvkkS5cupWvXruXuZ7VaiYqK4qGHHuKRRx6hSZMmjqavMZhVq1Yxa9YslixZwgMPPODUsr2ypgkL\nC+Pee+/l7NmzN6y3JTKr1crUqVN57rnnKhQMFI8IXL58OQ0aNGDkyJGcOnXKqXlr3Effvn15++23\nGTRoED//7HD3RrvxaC/nqKgo1q1bR1BQEL/++is7d+7kxx9/JDw8nMjISMxmM/fccw9ZWVn4+fnx\n+uuvV6rcTp060alTJ1q0aMGYMWOYM2eOrnG8lI4dOzJmzBj69OnD5s2badq0qctjeuzlGUBOTg4v\nv/wyx44d4+zZs0ycOJFHH32UvLw8MjIyCAgI4Pvvv6ewsJCRI0c6NNLvk08+YeLEiYwaNYq+ffve\nMH2Hns7De/jqq69ISUnh448/JikpqcrleWXrmTvZsmULEyZM4OzZsyQmJqKUYs+ePezYsYPo6Ghi\nY2OpVatWpcuzWq2li8VisfnXarWilMLPzw8/Pz9EBH9//9LXZddD8WXp9eX68WVjmc1mioqKSmvh\nyMhI6tatS+3atQkPD6dmzZrUrl2bhg0buuR/6Ans2bOHqVOn0rdvX+bMmVOlqSW1aCrJ+vXr+fe/\n/42fnx9du3alR48epKens2/fPrKysipdznUB+Pv7ExAQUPr65kVEfiMuW69F5AYh3SwsPz8/goOD\nCQwMJCvY5VDGAAAHwUlEQVQriwsXLnD+/HnS09PJyMggMzOTzMxMjh07xvTp0336QfG1a9f4y1/+\nwr59+5g3b57DDQRaNBoAfvjhB5544gkWLlzosaYVzuL7779nypQpXLt2zaFLbC0aTSkvvfQSP/30\nE7Nnz/ZYV35HsVgsnDhxggMHDrBx40by8vLYsWOHQ2W5RDQiMgCYCrQCOpXxCAgBPgXiKG6d+4dS\nakY5ZWjRuBmLxcLjjz/OxYsXeffddz3W8eVWKKW4dOkSBw4c4ODBg6VLVFRU6cDFQYMGOVyjuko0\nrQArxWaBZY01hgN9lVKDRSQUOAj0UEqdtlGGFo0BFBUVMXz4cLZs2cITTzxBr169qF27ttFpVUhR\nURFHjhxh79697N+/n71792K1WunUqRMJCQkkJCTQqVMn6tSp45R4FYnGFWaB54EaIuIP1ADMQLaj\ncTTOJyAggAULFrB27Vrmzp3L+++/T0xMDOHh4SilqFOnDjExMcTFxREfH09ERATXrl0jLy+P8PDw\n0vFNJ06cYM+ePcTExBAfH+/UDpTXrl3jwIED7Nmzh/3793PgwAEaN25Mt27dGDZsGL/73e9o0qSJ\nIR1xnW7hVLLuM4qdNsOAsUqpeeUcq2saD8BkMrFv3z5yc3MRES5cuMDBgwfZvHkzO3bsoG7dumRk\nZBAREUFOTg6NGzcmLy+PoqIi7r//fvbv38+5c+do1KgRkZGR9OvXj/j4eNLT05k3bx67d++moKCA\niIgIWrRoQUJCAklJSURERJTmkJ+fz549e9i2bRu7du3i5MmTxMfH0717d+677z7uvfdep9UilcHh\nmsYRs0ARGQqEArcDkcAmEVmrlDpha//yzAI17iM0NJQuXbrY3GY2mzl69CixsbEEBweTl5fHwYMH\nCQ0N5a677sLPr7gn1smTJzl//jyHDx9m3rx5TJkyhZCQEN544w2++OILwsPDyczMJDU1lYULFzJg\nwACeeeYZGjduzHfffcePP/5Iu3bt6Nu3LyNHjqRjx45ubagw2ixwLrBFKfVZyftPgJVKqS9tHKtr\nmmrKnj17+NOf/lTq2T106FCPuq9yaZNziWjGK6V2lbwfDbRXSv1BRGoA24GBSqkDNo7VotF4JC7p\n5VyeWSDFrWlBIrKfYsGk2BKMRuOt6IebGo0NvHI8jUbjqWjRaDR2okWj0diJFo1GYydaNBqNnWjR\naDR2okWj0diJFo1GYydaNBqNnWjRaDR2okWj0diJFo1GYydaNBqNnWjRaDR2okWj0dhJVQahJYvI\nIRHZKyLfiEitMtsmicgxETksIn2ck6pG4xlUpaZZDcQppe4GjgKTAESkNTAQaA08CMwVEafUaJU1\nPqjuMXzhM3hyDIe/zEqpNUopa8nbbUB0yeskYJFSqlApdRI4DnR2NE5ZPPWf6GkxfOEzeHIMZ93T\n/AH4d8nrhkDZ6cvOAnpWWI3PUGXfMxH5E2BWSi2soChtBKDxHcpOFmTvAgwHfgRCyqybCEws834l\n0KWc45Ve9OKpS3nf+6oYoD8IzKLY3DyjzPrWwEKK72MaAd8DLbTtjMZXqMqEkh8AQcCaEhPqn5RS\nzyulDorIEopnCygCnteC0fgShvqeaTTeiMf2CBCRFBG5WOLUeX1duQ9UnVF+mW3jRMQqIpGOll9R\nDBEZVfI5DojIu86OISKdRWS7iKSKyA4R6VTFGDEi8oOI/FyS8+iS9ZEiskZEjorIahFxyIy5gvKd\neb5txiizvfLnvCoNAa5cgG5APLC/zLregF/J6xnADGeWX7I+huLGixNApAs+Qy9gDRBY8v42F8RY\nT/HEWgCJwA9VjBFFsT83QDhwBLgLmAlMKFn/qqPno4LynXm+bcZw5Jx7bE2jlNoEXL1pXXkPVJ1S\nfgl/ASY4Wm4lYjwHvKOUKizZ57ILYpwHrv8q1wbOVTHGBaXUnpLX14BDFDfy/D/g7yW7/R3o78Ty\nGzr5fNuMUbLZrnPusaKpBGUfqDoFEUkCziql9jmz3JuIBbqLyFYRWS8iHV0QYyIwS0ROA8mUdHFy\nBiLShOKabRvQQCl1sWTTRaCBk8svi9POd9kYjpzzqrSeGUYlH6jaW2YYMJniS4LS1c4qvwwBQB2l\nVELJvcYSoJmTY3wCjFZKLS2ZUDiFGz+XQ4hIOPA1MEYplVN26j6llBKRKrUqlZT/VUn518qsd9r5\nLhuD4jlj7T7nXlfTlEyE+xAwxMlFNweaAHtF5ATFlwK7RKS+k+OcBb4BUErtAKwiUtfJMTorpZaW\nvP4KJ/T9E5FAigWzQCn1bcnqiyISVbL9duCSE8r/rEz5Tj3fNmI4dM69SjQlD1RfAZKUUvnOLFsp\ntV8p1UAp1VQp1ZTiL3cHpZTDX4Ry+Ba4H0BEWgJBSqkrTo5xXER6lLy+n+Je6A4jxVXKJ8BBpdSc\nMpv+CQwreT2M4s/mtPKdeb5txXD4nFelVcWVC7AISKd4dugzFF/THgNOAakly1wnlF9QUv7TN21P\no+qtZ7+JAQQCC4D9wC6gp5P/T08DHSm+J9gD/ATEVzHGfRRfyuwp879/kOI5Vb+nWJSrgdpOLD/R\nyefbZgxHzrl+uKnR2IlXXZ5pNJ6AFo1GYydaNBqNnWjRaDR2okWj0diJFo1GYydaNBqNnWjRaDR2\n8v8Bxt1uVMU3QdUAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.path as mpath\n",
"import matplotlib.patches as mpatches\n",
"\n",
"\n",
"# Extract first layer of features from shapefile using OGR\n",
"ds = ogr.Open('world.shp')\n",
"lyr = ds.GetLayer(0)\n",
"\n",
"\n",
"# Prepare figure\n",
"plt.ioff()\n",
"plt.subplot(1,1,1)\n",
"ax = plt.gca()\n",
"\n",
"\n",
"paths = []\n",
"lyr.ResetReading()\n",
"\n",
"lyr.SetAttributeFilter ( \" NAME = 'ANGOLA' \")\n",
"ax.set_xlim(11, 24.5 )\n",
"ax.set_ylim(-20, -2)\n",
"# Read all features in layer and store as paths\n",
"\n",
"for feat in lyr:\n",
"\n",
" for geom in feat.GetGeometryRef():\n",
" envelope = np.array( geom.GetEnvelope() )\n",
" # check if geom is polygon\n",
" if geom.GetGeometryType() == ogr.wkbPolygon:\n",
" codes = []\n",
" all_x = []\n",
" all_y = []\n",
" for i in range(geom.GetGeometryCount()):\n",
" # Read ring geometry and create path\n",
" r = geom.GetGeometryRef(i)\n",
" x = [r.GetX(j) for j in range(r.GetPointCount())]\n",
" y = [r.GetY(j) for j in range(r.GetPointCount())]\n",
" # skip boundary between individual rings\n",
" codes += [mpath.Path.MOVETO] + \\\n",
" (len(x)-1)*[mpath.Path.LINETO]\n",
" all_x += x\n",
" all_y += y\n",
" path = mpath.Path(np.column_stack((all_x,all_y)), codes)\n",
" paths.append(path)\n",
" # Add paths as patches to axes\n",
" for path in paths:\n",
" patch = mpatches.PathPatch(path, \\\n",
" facecolor='0.8', edgecolor='black')\n",
" ax.add_patch(patch)\n",
"\n",
"\n",
"\n",
"ax.set_aspect(1.0)\n",
"plt.show()"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 0
}