Vector Layers

1. Configure the library
2. Supported formats
3. Open datasources
4. Inspecting layers
5. Get layers features
6. Get layers fields
7. Get features geometries
8. Layer Spatial Referece Systems
9. Visualize the non geometric data in a deedle frame
10. Getting and Querying Open Street Map Data

References

The OGR part of GDAL/OGR Library provides a single vector abstract data model to the calling application for a variety of different supported formats including among others the well known ESRI Shapefile, but also RDBMSes, directories full of files, or even remote web services depending on the driver being used.

Configure the library

Initially it is necessary to configure the library and register all the format drivers that are desired. The FSharp.Gdal library accomplishes this calling the function Configuration.Init():

1: 

Configuration.Init()

Supported formats

The list of the supported formats is very big: to check that the library is correctly configured and all OGR drivers are properly registered we call:

1: 

Configuration.printOgrDrivers()

OGR 0: ESRI Shapefile
OGR 1: MapInfo File
OGR 2: UK .NTF
OGR 3: SDTS
OGR 4: TIGER
OGR 5: S57
OGR 6: DGN
OGR 7: VRT
OGR 8: REC
OGR 9: Memory
OGR 10: BNA
OGR 11: CSV
OGR 12: NAS
OGR 13: GML
OGR 14: GPX
OGR 15: LIBKML
OGR 16: KML
OGR 17: GeoJSON
OGR 18: Interlis 1
OGR 19: Interlis 2
OGR 20: GMT
OGR 21: GPKG
OGR 22: SQLite
OGR 23: ODBC
OGR 24: WAsP
OGR 25: PGeo
OGR 26: MSSQLSpatial
OGR 27: PostgreSQL
OGR 28: MySQL
OGR 29: PCIDSK
OGR 30: OpenFileGDB
OGR 31: XPlane
OGR 32: AVCBin
OGR 33: AVCE00
OGR 34: DXF
OGR 35: Geoconcept
OGR 36: GeoRSS
OGR 37: GPSTrackMaker
OGR 38: VFK
OGR 39: PGDump
OGR 40: OSM
OGR 41: GPSBabel
OGR 42: SUA
OGR 43: OpenAir
OGR 44: PDS
OGR 45: WFS
OGR 46: HTF
OGR 47: AeronavFAA
OGR 48: Geomedia
OGR 49: EDIGEO
OGR 50: GFT
OGR 51: GME
OGR 52: SVG
OGR 53: CouchDB
OGR 54: Idrisi
OGR 55: ARCGEN
OGR 56: SEGUKOOA
OGR 57: SEGY
OGR 58: ODS
OGR 59: XLSX
OGR 60: ElasticSearch
OGR 61: PDF
OGR 62: Walk
OGR 63: CartoDB
OGR 64: SXF

Open datasources

Next we need to open the input OGR datasource. Despite the varietry of formats the datasource name is always a single string.

In this case we will open a shapefile provided by the Italian National Institute of Statistics (Istat) containing the administrative limits of italian regions.

1: 
2: 
3: 
4: 
5: 

open OSGeo

let italyShp = __SOURCE_DIRECTORY__ + @"\data\reg2001_g\reg2001_g.shp"
let driver = OGR.Ogr.GetDriverByName("ESRI Shapefile")
let ds = driver.Open(italyShp, 0)

OSGeo.OGR.DataSource

An OGR.DataSource contains a set of OGR.Layers each of which in turn contains a set of OGR.Features and an OGR.FeatureDefn.

Globally we can though of the datasource as a database, a layer as a table in the database and the layer's features as its records while the feature definition provides its columns.

From the output above this structure is not visbile. The aim of FSharp.Gdal is to improve the use of the standard OGR/GDAL library within the fsharp interactive shell and to provide utility functions to make it easier (and more idiomatic) writing F# scripts to process geospatial data.

To this aim the first thing that we want is to get immediately from the output an overview of the objects we're dealing with and an hint to successive inspections. The first of this functions is datasourceInfo implemented in the FSharp.Gdal.Vector module (below we will use more of these functions the implementation of which can be viewed directly in the source code of the Vector module ).

1: 

let dsInfo = ds |> Vector.datasourceInfo

[{Index = 0;
  Layer = {Features = 20;
           Geometry = wkbPolygon;
           Fields = ["COD_REG"; "REGIONE"; "POP2001"];};}]

Not only we can call directly the datasourceInfo function but we can also use it to define a custom printer for the OGR.DataSource:

1: 

fsi.AddPrintTransformer (fun (datasource:OGR.DataSource) -> box (datasource |> Vector.datasourceInfo))

In the sequel we will inspect the layer's structure and define a custom printer also for the OGR.Layer after that whenever we will open again a new datasource we will get all the infos describing its structure.

Inspecting layers

As pointed out above an OGR.DataSource can potentially have many layers associated with it. The number of available layers can be queried with GetLayerCount() and individual layers fetched by index using GetLayerByIndex.

1: 
2: 
3: 

let layersCount = ds.GetLayerCount()

printfn "There is %i layer" layersCount

There is 1 layer

1: 

let layer0 = ds.GetLayerByIndex(0)

instead of using these low level GDAL/OGR methods however the Vector module provides the function layers which returns the layers in a Vector DataSource as an F# list.

1: 

let layers = ds |> Vector.layers

[OSGeo.OGR.Layer]

As we can see the output is still not very informative so we call the function layerInfo to immediately inspect the layer's content:

1: 
2: 
3: 

let layerInfo = 
    layers
    |> List.map (fun ly -> ly |> Vector.layerInfo)

1: 

printfn "%s" (layerInfo.ToString())

[{Features = 20;
Geometry = wkbPolygon
Fields = ["COD_REG"; "REGIONE"; "POP2001"]}]

This is better: now we can see that the layer has 20 features (records) in it, is a layer of geometry type polygon and all features in it provide three non geometric informations (attributes or fields): COD_REG, REGIONE and POP2001.

To complete the configuration of the F# interactive shell we use the function to define the custom printer for the OGR.Layer type:

1: 

fsi.AddPrintTransformer (fun (layer:OGR.Layer) -> box (layer |> Vector.layerInfo))

Get layers features

A layer contains a set of features that can be viewed as records of a database table. To get them just call features to obtain an F# list:

1: 
2: 

let layer = layers |> List.head
let features = layer |> Vector.features

Get layers fields

The feature definition is globally associated to the containing layer. Here the function to use is fields which returns a tuple of all the index, name and type of the fields in the layer's definition.

1: 
2: 
3: 
4: 

let attributes = layer |> Vector.fields

for fieldIndex, fieldName, fieldType in attributes do
    printfn "Field Index = %i; Field Name = %s; Field Type = %A" fieldIndex fieldName fieldType

Field Index = 0; Field Name = COD_REG; Field Type = OFTInteger
Field Index = 1; Field Name = REGIONE; Field Type = OFTString
Field Index = 2; Field Name = POP2001; Field Type = OFTInteger

In this case the attribute "Regione" is the region's name so we can iterate over all features to get all the italian regions' names quering the field with index 1:

1: 
2: 
3: 

for feat in features do 
    let regionName = feat.GetFieldAsString(1)
    printfn "%s" regionName

PIEMONTE
VALLE D'AOSTA
LOMBARDIA
TRENTINO-ALTO ADIGE
VENETO
FRIULI VENEZIA GIULIA
LIGURIA
EMILIA-ROMAGNA
TOSCANA
UMBRIA
MARCHE
LAZIO
ABRUZZO
MOLISE
CAMPANIA
PUGLIA
BASILICATA
CALABRIA
SICILIA
SARDEGNA

Get features geometries

In addition to attributes a feature contains a geometry field. In this example each feature corresponds to an italian region. To graphically plot them we can collect all features' geometries in a geometry collection and call the Plot utility (see Plot Geometry) to get a map of Italy divided by regions:

1: 
2: 
3: 
4: 
5: 
6: 

let geomcol = new OGR.Geometry(OGR.wkbGeometryType.wkbGeometryCollection)

for feat in features do 
    geomcol.AddGeometry(feat.GetGeometryRef()) |> ignore

let italyRegionsPlot = geomcol |> Plot

Spatial reference system

A geometry is always defined in a coordinate system but in the case of geospatial data we talk more properly of Spatial Referece Systems:

1: 

let _,sr = layer.GetSpatialRef().ExportToWkt()

1: 

printfn "%s" sr

PROJCS["ED_1950_UTM_Zone_32N",GEOGCS["GCS_European_1950",DATUM["European_Datum_1950",SPHEROID["International_1924",6378388.0,297.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

Visualize the non geometric data in a deedle frame

Deedle is an F# library for data exploration and manipulation. FSharp.Gdal wants to be deedle friendly and while it does not reference directly this library it provides the function toValues to immediately convert vector geospatial data in a deedle frame:

1: 
2: 
3: 
4: 
5: 
6: 

open Deedle

let fmItalyRegion = 
    layer
    |> Vector.toValues
    |> Frame.ofValues

As we can see the data converted in a deedle frame are more readable and we can take advantage of the library capabilities to make aggregations, pivoting and scientific and statistical calculations on the data.

Getting and Querying Open Street Map Data

As pointed in the beginning GDAL/OGR library has the capabilites of access not only file datasources as shapefiles, but also RDBMSes, web services, etc. In this section we can see how to use it to get data from the OpenStreetMap.

OpenStreetMap is one of largest (probably the largest accessible one) source of geospatial data. It provides not only a free map accessible from the website Openstreetmap.org but also an Editing API for fetching and saving raw geodata from/to the OpenStreetMap database and an Overpass API which provides read-only API access.

In this section we will see how to

1. fetch data from the Overpass API using the methods in the FSharp.Data library
2. save the data to a file (standard .NET methods)
3. querying the saved data mixing the GDAL/OGR library's methods with F# and in particular with the Deedle library

The first thing is to download the data. In this case we want all the streets in Malesco (one of the "biggest" city centre at the border of the Val Grande National Park).

I won't go much into the deatils of the query syntax here which can be found on the wiki site Overpass API but just to be clear we are querying for all the ways (streets, waterway, etc.): way[name=*];in the bounding box containing Malesco city: [bbox=8.4872,46.1220,8.5140,46.1331]

1: 
2: 
3: 
4: 
5: 
6: 
7: 
8: 

open FSharp.Data

//bbox = left, bottom, right, top (min long, min lat, max long, max lat)

let baseUrl = "http://overpass.osm.rambler.ru/cgi/xapi?way[name=*][bbox=8.4872,46.1220,8.5140,46.1331]"

// Download the content
let osm = Http.RequestString(baseUrl)

Then we save the downloaded data to a new file:

1: 
2: 
3: 
4: 
5: 
6: 
7: 

let fileName = __SOURCE_DIRECTORY__ + @"\data\malescoWays.osm"

let save  (fileName:string) (data:string) = 
    use w = new System.IO.StreamWriter(fileName)
    w.Write(data)

osm |> save fileName

Now we can inspect the data using the GDAL/OGR library.

First we open the datasource and get the layers inside it:

1: 

let osmDataSource = OGR.Ogr.Open(fileName, 0)

"[{Index = 0;
  Layer =
   {Features = 9;
    Geometry = wkbPoint;
    Fields =
     ["OSM_ID"; "NAME"; "BARRIER"; "HIGHWAY"; "REF"; "ADDRESS"; "IS_IN"; "PLACE";
      "MAN_MADE"; "OTHER_TAGS"];};};
 {Index = 1;
  Layer =
   {Features = 55;
    Geometry = wkbLineString;
    Fields =
     ["OSM_ID"; "NAME"; "HIGHWAY"; "WATERWAY"; "AERIALWAY"; "BARRIER";
      "MAN_MADE"; "OTHER_TAGS"];};};
 {Index = 2;
  Layer = {Features = 0;
           Geometry = wkbMultiLineString;
           Fields = ["OSM_ID"; "NAME"; "TYPE"; "OTHER_TAGS"];};};
 {Index = 3;
  Layer =
   {Features = 3;
    Geometry = wkbMultiPolygon;
    Fields =
     ["OSM_ID"; "OSM_WAY_ID"; "NAME"; "TYPE"; "AEROWAY"; "AMENITY";
      "ADMIN_LEVEL"; "BARRIER"; "BOUNDARY"; "BUILDING"; "CRAFT"; "GEOLOGICAL";
      "HISTORIC"; "LAND_AREA"; "LANDUSE"; "LEISURE"; "MAN_MADE"; "MILITARY";
      "NATURAL"; "OFFICE"; "PLACE"; "SHOP"; "SPORT"; "TOURISM"; "OTHER_TAGS"];};};
 {Index = 4;
  Layer = {Features = 0;
           Geometry = wkbGeometryCollection;
           Fields = ["OSM_ID"; "NAME"; "TYPE"; "OTHER_TAGS"];};}]"

We have five layers respectively of wkbPoint, wkbLineString, wkbMultiLineString, wkbMultiPolygon, wkbGeometryCollection geometry type and as we can see the third and fifth layers are empty.

We are interested in the streets so let's convert the second wkbLineString layer in a deedle frame to see its content:

1: 
2: 
3: 
4: 
5: 
6: 
7: 
8: 
9: 

let osmLayers = 
    osmDataSource 
    |> Vector.layers

let fmMalescoWays = 
    osmLayers
    |> Seq.item 1
    |> Vector.toValues
    |> Frame.ofValues

We can see also that this layer contains ways of different kind: only the records with the a value in the column HIGHWAY are streets so let's see the list of its distinct values:

1: 
2: 
3: 
4: 

let highways = 
    fmMalescoWays.GetColumn<string>("HIGHWAY").Values 
    |> Seq.distinct
    |> List.ofSeq

["secondary"; "tertiary"; "unclassified"; "residential"; ""; "service"; "track";
 "path"; "footway"]

We can now filter just the streets of residential type:

1: 
2: 
3: 
4: 
5: 
6: 

let fmMalescoStreets = 
    fmMalescoWays 
    |> Frame.filterRowValues 
        (
            fun row -> row.GetAs<string>("HIGHWAY") = "residential"
        )

References

• Deedle
• Italian National Institute of Statistics (Istat)
• OGR API Tutorial
• Openstreetmap.org
• Overpass API
• Python GDAL/OGR Cookbook

• Fsharp.Gdal
• Home page
• Get Library via NuGet
• Source Code on GitHub
• License
• Release Notes
• Getting started
• Sample tutorial
• Documentation
• API Reference