The problem...

If you publish a SQL Server (in this case 2014) table or view with a GUID field using GeoServer (version 2.10), it will fail to make  a valid WFS because the &request=DescribeFeatureType response will not include these layers.

This manifests itself in various ways depending on the client, but in QGIS on trying to load the layer into the map we get this in the "Log Messages Panel":

Analysis of DescribeFeatureType response failed for url  srsname='EPSG:27700' typename=xxx:yyyurl='https://mysite.com/myworkspace/wfs' version='1.0.0' table="" sql=: it is probably a schema for Complex Features

which had me puzzled for quite some time!

The solution...

Changing the data source view to include this field with

CAST(myGUID AS varchar(36)) AS myGUID

then hitting "Reload feature type" in the GeoServer layer configuration screen fixes it. This changes the field data type from "UUID" to "string" in that screen. Or just omit the field entirely.

My guess is this will not help transactional WFS, but I have not tested this. Overall, GUIDs are still best avoided in GIS-land as they are so poorly supported by the database drivers.

R is a cross platform statistical package which is becoming extremely widely used.  It is modular and so there are all sorts of add ins available, including a number of sophisticated tools for spatial analysis some of which run considerably faster than Arc / MapInfo.  Organisations are also starting to move to this for spatial analysis as its cheaper but more relevantly doesn’t suffer from the endless version changes that ArcGIS is particular forces on people.

R is generally used for a wide range of statistics – but also has the option to generate publication quality graphs (pdfs or images)

Getting started

R can be downloaded here: http://www.r-project.org/

Its also suggested (on Windows) to use r studio - http://www.rstudio.com/ which gives a more friendly IDE

R is a scripted tool – so it can be a steep learning curve but there are endless examples on the web to help get started.

Spatial example using OGR

This syntax is a little obscure – but potentially of considerable use for a wide range of spatial data.  Reading from shapefiles (or any other OGR supported datasource is easy / possible)

Read data from SQL server

#download and load rgdal libraries
install.packages("rgdal")
library(rgdal)

#Use ogrinfo to get information about a shape file called stations
ogrInfo("C:/ToFile/RNLI/data", "Stations")

#List the spatial layers available in a SQL server database
#Note you must have the standard OGR tables geometry_columns and spatial_ref_sys present and populated in your database
#Some people say you need a DSN as well – not sure you do …
ogrListLayers("MSSQL:server=SQL2008;database=CMSI-NE;trusted_connection=yes")

#Read a table of polygon objects into a new sitesPoly object
#Then plot them on a map
#N.B.  You can’t have binary data in your table – so if you do create a view or somesuch
sitesPoly <- readOGR(dsn='MSSQL:server=SQL2008;database=CMSI-NE;trusted_connection=true', layer='SitesGeom')
plot(sitesPoly)

#Select just one record from the layer
oneSite = subset(sitesPoly, sitesPoly$MI_PRINX == 420)

#Make sure the layer is assigned the correct projection
proj4string(oneSite) = CRS("+init=epsg:27700")

#Transform the layer to another projection
oneSite_latlong = spTransform(oneSite, CRS("+init=epsg:4326"))

#Load another library which contains some world boundaries
#Load up coast and countries outlines
library(rworldmap)
data(coastsCoarse)
data(countriesLow)

#plot countries outline limiting extent to Europe
plot(countriesLow, xlim = c(-10, 10), ylim = c(45, 60), legend = F)

#Add sites onto existing map setting colours

proj4string(sitesPoly_ll) = CRS("+init=epsg:27700")
sitesPoly_ll = spTransform(sitesPoly, CRS("+init=epsg:4326"))
plot(sitesPoly_ll, col = "red", border="red", Add=T)

Use a different library to do some plotting

install.packages("maps")
install.packages("mapdata")
library(maps)
library(mapsdata)

#Select polygons making up UK and plot
#Then plot a point on top
map('worldHires',
    c('UK', 'Ireland', 'Isle of Man','Isle of Wight','Wales:Anglesey'),
    xlim=c(-11,3), ylim=c(49,60.9))
points(-1.615672,54.977768,col=2,pch=18)

In another couple of areas, the WMS requests to MapServer simply failed with this error sequence:

[Sat Jan 21 09:30:02 2012].785000 msDrawRasterLayerLow(OS50KColourRaster): entering.
[Sat Jan 21 09:30:02 2012].791000 msResampleGDALToMap in effect: cellsize = 5.000000
[Sat Jan 21 09:30:02 2012].792000 msDrawGDAL(OS50KColourRaster): using RAW_WINDOW=3509 0 491 652, dst=0,0,491,652
[Sat Jan 21 09:30:02 2012].792000 msDrawGDAL(): red,green,blue,alpha bands = 1,0,0,0
[Sat Jan 21 09:30:02 2012].891000 msResampleGDALToMap in effect: cellsize = 5.000000
[Sat Jan 21 09:30:02 2012].892000 msDrawGDAL(OS50KColourRaster): using RAW_WINDOW=3509 3348 491 652, dst=0,0,491,652
[Sat Jan 21 09:30:02 2012].892000 msDrawGDAL(): red,green,blue,alpha bands = 1,0,0,0
[Sat Jan 21 09:30:02 2012].998000 msResampleGDALToMap in effect: cellsize = 5.000000
[Sat Jan 21 09:30:02 2012].998000 msDrawGDAL(OS50KColourRaster): using RAW_WINDOW=0 0 812 652, dst=0,0,812,652
[Sat Jan 21 09:30:02 2012].998000 msDrawGDAL(): red,green,blue,alpha bands = 1,0,0,0
[Sat Jan 21 09:30:03 2012].122000 msResampleGDALToMap in effect: cellsize = 5.000000
[Sat Jan 21 09:30:03 2012].122000 msDrawGDAL(OS50KColourRaster): using RAW_WINDOW=0 3348 812 652, dst=0,0,812,652
[Sat Jan 21 09:30:03 2012].122000 msDrawGDAL(): red,green,blue,alpha bands = 1,0,0,0
[Sat Jan 21 09:30:03 2012].148000 GetBlockRef failed at X block offset 0, Y block offset 499: Unable to access file. GDALDatasetRasterIO() failed: drawGDAL()
[Sat Jan 21 09:30:03 2012].149000 msDrawMap(): Image handling error. Failed to draw layer named 'OS50KColourRaster'.
[Sat Jan 21 09:30:03 2012].149000 freeLayer(): freeing layer at 024265B0.

again the underlying images were apparently fine.

And of course this has left gaps in the tile cache.

So, in order to overcome these issues I am setting up the WMS using GeoServer to see if that can fill the gaps.

My starting point is a folder containing the TIFF images, each with a .TFW world file.

Gotchas:

There are instructions for setting up a raster WMS in GeoServer here, but here is my summary with details of where things can go wrong:

• Create a workspace, in my case “ESDM_UK_BaseMaps” with a namespace of “https://www.ordnancesurvey.co.uk/”
• Create a store – this will equate with one type of map, in this case the Landranger 1:50,000 colour raster.
• Choose ImageMosaic as the type of store
• Select the workspace, and give the store a name, in this case “OS50K”
• For the connection parameters, enter the folder that contains the images, in my case “file:D:\MapData\OrdnanceSurvey\50KRasterColour”
• GeoServer will then automatically scan the folder and create a shapefile index for the tiles (or “granules” as they bizarrely call them). At this point I was confronted with this error message “Could not list layers for this store, an error occurred retrieving them: Unable to acquire a reader for this coverage with format: ImageMosaic”.

This is a commonly observed error by all accounts with lots of fairly grumpy exchanges among GeoServer users and developers about how to reproduce it and what might be the causes. It seems it can be raised by a raft of different issues. So I returned to a close reading of the instructions and also found Adrian Walker’s useful blog post about configuring OS maps in GeoServer.

It turns out I had two problems:

First the .TFW extension must be lower case. So in a command prompt run

ren *.TFW *.tfw

(note that the image file extension does not have to be lower case)

Second, GeoServer needs a .prj file for every image (not required by MapServer).

These can be generated using a python script (thanks to Adrian Walker again for this tip):

import glob
content = 'PROJCS["OSGB 1936 / British National Grid", GEOGCS["OSGB 1936", DATUM["OSGB_1936", SPHEROID["Airy 1830",6377563.396,299.3249646, AUTHORITY["EPSG","7001"]], AUTHORITY["EPSG","6277"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.01745329251994328, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4277"]], UNIT["metre",1, AUTHORITY["EPSG","9001"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",49], PARAMETER["central_meridian",-2], PARAMETER["scale_factor",0.9996012717], PARAMETER["false_easting",400000], PARAMETER["false_northing",-100000], AUTHORITY["EPSG","27700"], AXIS["Easting",EAST], AXIS["Northing",NORTH]]'
tifs = glob.glob('*.tif') 
for tif in tifs:  
    prj = tif.split('.')[0] + '.prj'  
    file = open(prj,'w')  
    file.writelines(content)  
    file.close()

After resolving these issues, it was necessary to return to step 2 above to create the store, but I got the error again! OK, GeoServer had created the index shapefile already, and it was bad due to the missing prj and tfw files; unfortunately GeoServer will not overwrite this file so I had to manually delete the files (the shapefile is named after the store, and comprises seven files with various extensions).  After that, creating the store worked fine.

Hint: if you examine the .prj file in a text editor and it looks like this, then you need to bin it and start again:

LOCAL_CS["Generic cartesian 2D",
LOCAL_DATUM["Unknow", 0],
UNIT["m", 1.0],
AXIS["x", EAST],
AXIS["y", NORTH]]

So on clicking “Save” for the new store, it showed a twirly for a couple of minutes while creating the index shapefile, and then transfers you to the New Layer page, showing the new layer that can be published (it names the layer after the folder containing the images). From there it is a simple matter to complete the configuration of the layer.

Note: I think I could have configured a single store one folder higher up, and this would have included each subfolder of TIFFs as a separate layer.

Outcome

Once I had finally got my GeoServer WMS up and running it was time to see whether could fills the gaps in my tilecache…. and I’m very pleased to say yes it could!

I’ve posted separately about performance and image quality from the two servers.

1. setting up and accessing SQL Azure databases
2. using SQL Azure with ArcGIS + ArcSquirrel and MapLink
3. using SQL Azure with Quantum GIS and MapServer.

Azure Management Portal

Although I’ve never set up any cloud services before (other than Arvixe web hosting) it was incredibly easy to setup a trial SQL Azure subscription (“Windows Azure Platform Introductory Special”), create a database, and add some firewall rules to control access. I didn’t dig deep, but the “Management Portal” seems nice and simple:

It helpfully tells you your connection string for ADO.Net, ODBC and PHP (passwords not included!):

See this page for details of the Windows Azure Platform Introductory Special including what you get for no charge and the pricing for the rest.

MapLink, ArcSquirrel and Azure

Hartmut connected up our MapLink to the database and pushed some spatial data into my Azure database from ArcGIS – 20,000 polygons took about 75 minutes.  Hartmut had the layers drawing in ArcMAP, performing as well as in a normal SQL Server instance apparently, with only minor mods to our code to make the initial connection. However editing will not work without some more tweaks - MapLink and ArcSquirrel rely on some SQL metadata that is a little different on Azure.

SQL Server Management Studio and MS Access

Connecting with SSMS is simple, however you get reduced functionality – for example you cannot right-click > edit top 200 records or anything, so quick edits of a small table need to be done with SQL or in MS Access or some other client. Fortunately connecting with MS Access (ODBC + SQL Native Client 10) is also simple, and  the Azure database behaves just like a local SQL instance – remarkable querying speed (all things considered), though slower inserts (10,000 records took something like 15 minutes). NB you cannot directly link a spatial table into Access because it doesn’t recognize the geometry data type; but to update attributes you can link a view.

Also, in SSMS you get no designer tools, e.g. for creating a view; everything has to be done with T-SQL.

SQL Azure Database Manager

After writing the above I found the database manager in the Azure control panel, which has functionality similar to SSMS, including the ability to edit the tables*:

* spatial data types not supported for view or edit.

You get a designer for creating/modifying tables, but views and stored procedures are T-SQL only.

MapServer, OGR and Azure

Next I wanted to try serving the data through MapServer, and viewing/editing the data in Quantum GIS or MapInfo.  MapServer and QGIS require the “geometry_columns” and “spatial_ref_sys” tables, and the easiest way to create these is to push some data in with ogr2ogr.

So I tried this from one of our London servers. I loaded a few hundred linear features, but got four errors like “INSERT command for new feature failed… incorrect syntax…”. On checking the results, a) while the geometries had loaded fine, all the attribute fields were NULL, and b) the spatial_ref_sys table was empty.  I had four attribute fields, so my guess is that each syntax error was one field update across all records using a join on the incoming data.

Next I hooked up a MapServer WMS to my SQL Azure tables. MapServer had no trouble connecting to my Azure tables, but initially wouldn’t draw anything…  on investigating I found that ogr2ogr had not correctly set the SRID of the geometries, so I had to update them “manually” with:

update esdm_waters set ogr_geometry.STSrid = 4326

(which, incidentally, took 1 minute 52 seconds for 788 rows, demonstrating that updating an Azure database ain’t quick).

Clearly OGR needs a little work to do on SQL Azure support, but this was enough to get me going.

The WMS worked, so to challenge it a little more I set up another WMS onto a table with about 8000 sub-tidal habitat polygons, showing red in this OpenLayers map:

And here showing the rivers WMS as well (blue lines):

Performance of these two layers was perfectly acceptable, though actually about 100% slower than the same WMS working from a local SQL Server 2008 Express instance. It would be interesting to see what it could do with MapServer also running on Azure.

I also tried a WFS – again, performance was acceptable though significantly slower than running with a local SQL instance (though I didn’t do strict comparisons).

Quantum GIS and SQL Azure

I found this really helpful blog post to get started: http://woostuff.wordpress.com/2011/03/13/opening-ms-sql-spatial-in-qgis/

Adding the layers requires some python script, but I’m sure someone will write a GUI plugin for QGIS soon.

First make sure you are running a version of QGIS that has GDAL/OGR 1.8 libraries – I’m on 1.7.0 which is OK.

Open the Plugins menu > Python Console

and type in these two commands (modified for your data):

uri = "MSSQL:server={serverName};database={databaseName};tables={tableName};trusted_connection=yes"
qgis.utils.iface.addVectorLayer(uri,'{yourLayerNameHere}','ogr')

In my case:

uri = "MSSQL:server=tcp:a0awtn7gg9.database.windows.net;database=ArcSquirrel1;tables=esdm_waters;Uid=exegesis@a0awtn7gg9;Pwd=myPassword;"
qgis.utils.iface.addVectorLayer(uri,'Rivers','ogr')

Without any fuss my rivers dataset appeared:

And I can create and edit features and their attributes, here adding a test line off the Hebrides:

And to prove it, here’s my new feature in OpenLayers…

And finally, here we have the habitat polygons (EPSG:27700) colour coded by habitat type, as well as the rivers (EPSG:4326), in QGIS. Polygon editing also worked fine, once I’d worked out why my new polygons kept vanishing (doh, classified symbols).

Lots to think about, and I look forward to editing the same data with ArcGIS and ArcSquirrel soon…

...

Imports OSGeo.OGR

...

Dim drv As OSGeo.OGR.Driver = OSGeo.OGR.Ogr.GetDriverByName("PostgreSQL")
Dim connStr as string = "PG:dbname='myDB' host='192.168.etc.etc' port='5433' user='foo' password='bar'"
Dim ds As OSGeo.OGR.DataSource  = drv.Open(connStr, 0)

EXPRESSION ('[HAB_TYPE]' =~/^A1.1/) 

 If you try this with mapserver 6 however (especially on large datasets) you can get several mapserver threads starting - each of which ended up consuming 2GB RAM on the server - and never finishing

Changing the expression to:

EXPRESSION ('[HAB_TYPE]'  ~ '^A1.1') 

 fixes it.

This is a nasty issue as most changes to the map file are "breaking" changes so it just doesn't work.  This does run - but appears to kill the server.