1. Download the 64-bit version on the MS Access ODBC driver (AccessDatabaseEngine_x64.exe): http://www.microsoft.com/en-gb/download/details.aspx?id=13255
2. If you don’t have 32-bit MS Office installed you can just run the exe, however if you do have 32-bit office installed then you need to start an admin command prompt and run the installer in passive mode as follows:

AccessDatabaseEngine_x64.exe /passive

3. For QGIS 2: Then you need to edit the qgis-ltr.bat (normally in C:\Program Files\QGIS x.xx\bin, in previous versions used to be qgis.bat) and add the following lines (before the ‘start’ command)

set OGR_SKIP=ODBC
set PGEO_DRIVER_TEMPLATE=DRIVER=Microsoft Access Driver (*.mdb, *.accdb);DBQ=%%s

4. For QGIS 3: In QGIS Settings panel -> Options | System | Environment add the following two variables:

   variable name: PGEO_DRIVER_TEMPLATE 
   value: DRIVER=Microsoft Access Driver (*.mdb, *.accdb);DBQ=%s
   Or if being used on a shared server add a system environment variable 
   PGEO_DRIVER_TEMPLATE=DRIVER=Microsoft Access Driver (*.mdb, *.accdb);DBQ=%s

5. Now you should be able to start QGIS and simply drag the geodatabase mdb file onto the QGIS workspace and it will open.

Aside: I should say that what follows is not what we normally do in our real applications – we usually use ASP.Net .ashx handlers to relay requests through the root of our web sites, which also helps to avoid cross-domain problems with GetFeatureInfo responses and WFS, but that’s another story! This post is in reaction to Dennis’s request for help, and because I wanted to check this approach is viable.

I’ve posted before on how to use URL Rewrite to give GeoServer a friendly URL on IIS, and we’re going to use the same technique here, so first ensure you have the necessaries installed.

You should also already have a working WMS/WFS – we created this one in the previous post:

http://mapservertest.esdm.co.uk/scripts/mapserv.exe?map=D:\Websites\mapservertest\map\mapfiletest.map

which was fine, but a) the URL reveals some information about the structures on the web server, and b) the URL is a bit ugly.

So decide what address you want for your MapServer service (remembering this will be specific right down to the map file). I’m using “mapservertestwrapper.esdm.co.uk” which I have pointed at the web server in question.

Create a new web site in IIS with this URL as a host header (you’ll first need to create an empty “physical” folder for the site to point at – add a readme.config file explaining what it’s for).

Now open URL Rewrite and create a reverse proxy inbound rule, entering the real service address as the forwarding address:

OK now we should be able to see our MapServer on the new URL. However, when we visit http://mapservertestwrapper.esdm.co.uk/ we get this feedback in the browser:

msLoadMap(): Regular expression error. MS_DEFAULT_MAPFILE_PATTERN validation failed. msEvalRegex(): Regular expression error. String failed expression test.

The fix is to re-edit the inbound rule, and remove the final “/” in the “rewrite URL”, going from:

http://mapservertest.esdm.co.uk/scripts/mapserv.exe?map=D:\Websites\mapservertest\map\mapfiletest.map/{R:1}

to:

http://mapservertest.esdm.co.uk/scripts/mapserv.exe?map=D:\Websites\mapservertest\map\mapfiletest.map{R:1}

Now we should get sensible MapServer feedback in the browser (conditioned by a template specified in the map file).

Now let’s see if we can do a GetCapabilities request… turns out our URL has to be this:

http://mapservertestwrapper.esdm.co.uk/?service=wms&version=1.3.0&request=getcapabilities

Yes it works, but a) that’s a slightly odd URL with “/?”, and b) the response contains the original service addresses, again revealing information about our server.

To fix a), we go back and edit our inbound rule again, and change the “pattern” from “(.*)” to “map1(.*)” (without the quotes of course). Now our service URL becomes:

http://mapservertestwrapper.esdm.co.uk/map1?service=wms&version=1.3.0&request=getcapabilities

Ah yes that’s much nicer, and now of course we see we can manage lots of .map files within one web site simply by setting up in inbound rule for each, with an appropriate pattern and target.

To fix b) we can edit the metadata in the .map file to set this:

"wms_onlineresource"    "http://mapservertestwrapper.esdm.co.uk/map1"

This works fine, returning this within the GetCapabilities:

<OnlineResource xlink:href="http://mapservertestwrapper.esdm.co.uk/map1?" xmlns:xlink="http://www.w3.org/1999/xlink"/>

And I can confirm this creates a working map in QGIS, with map requests like this visible in Fiddler:

http://mapservertestwrapper.esdm.co.uk/map1?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=52801.44371538766426966,3347.3560407145996578,657116.82161524775438011,428145.51874091033823788&CRS=EPSG:27700&WIDTH=1020&HEIGHT=717&LAYERS=grid&STYLES=&FORMAT=image/png&DPI=96&MAP_RESOLUTION=96&FORMAT_OPTIONS=dpi:96&TRANSPARENT=TRUE

And a quick test shows that info clicking (GetFeatureInfo) and legends work too:

e.g. http://mapservertestwrapper.esdm.co.uk/map1?version=1.3.0&service=WMS&request=GetLegendGraphic&sld_version=1.1.0&layer=grid&format=image/png&STYLE=default

gives this:

There you go Dennis, we have working OGC services from MapServer with nice clean URLs, a simple architecture for adding new services with their own map files, and no leaking of details about our directory structures.

Footnote…

I’m not completely keen on the last bit, having to change the wms_onlineresource, as this means the original service is now broken (depending on whether the mapping client makes us of this, which QGIS likes to). Another way should be  to create an outbound URL Rewrite rule to rewrite any references to the original URL with the new one. However, after a bit of pfaffing with this, I’m pretty sure there’s a bug in URL Rewrite that means outbound rewriting doesn’t work for attributes containing colons, and ours is “xlink:href”. Incidentally setting up the outbound rules for custom tags is a bit of an adventure through the UI, and there is a bit of escaping to do for the backslashes and question mark in the URL (highlighted below) as URL Rewrite is doing regular expressions so you might as well get down and dirty in the web.config located in the “dummy” folder we created for the new web site we created earlier. So this is my outbound rule that doesn’t work, but should when someone fixes the bug!

<?xml version="1.0" encoding="UTF-8"?>
<configuration>
    <system.webServer>
        <rewrite>
            <rules>
                <rule name="ReverseProxyInboundRule2" stopProcessing="true">
                    <match url="map1(.*)" />
                    <action type="Rewrite" url="http://mapservertest.esdm.co.uk/scripts/mapserv.exe?map=D:\Websites\mapservertest\map\mapfiletest.map{R:1}" />
                </rule>
            </rules>
            <outboundRules>
                <rule name="ReverseProxyOutboundRule1" preCondition="ResponseIsXML" enabled="true">
                    <match filterByTags="CustomTags" customTags="OGC GetCapabilities" pattern="^http(s)?://mapservertest.esdm.co.uk/scripts/mapserv.exe\?map=D:\\Websites\\mapservertest\\map\\mapfiletest.map(.*)" />
                    <action type="Rewrite" value="http{R:1}://mapservertestwrapper.esdm.co.uk/map1{R:2}" />
                </rule>
                <preConditions>
                    <preCondition name="ResponseIsHtml1">
                        <add input="{RESPONSE_CONTENT_TYPE}" pattern="^text/html" />
                    </preCondition>
                    <preCondition name="ResponseIsXML">
                        <add input="{RESPONSE_CONTENT_TYPE}" pattern="^text/xml" />
                    </preCondition>
                </preConditions>
                <customTags>
                    <tags name="OGC GetCapabilities">
                        <tag name="OnlineResource" attribute="xlink:href" />       
                    </tags>
                </customTags>
            </outboundRules>
        </rewrite>
    </system.webServer>
</configuration>

The red bit is the bit affected by the bug I think, unless someone can point out how this should be escaped/expressed to make it work?

Reasons for this might include:

• Re-projecting the map
• Accessing a secured WMS/WFS and then protecting it a different way
• Serving a WMS / WFS map out in a different way

This is all relatively simple using standard mapserver .MAP file if the WMS / WFS services you are calling to are unsecured.  Then its simply a case of a layer like:

 
LAYER

NAME "MyLayer"
DEBUG 5
  TYPE RASTER
  STATUS ON
  CONNECTION http://somedomain.com/WMSServer
  CONNECTIONTYPE WMS
  METADATA
    "wms_title"              "OSI"
    "wms_srs"             "EPSG:27700 EPSG:4326 EPSG:54004 EPSG:3978 EPSG:3857 EPSG:900913"
    "wms_name"            "1,2,3,4,5,6,7,8,9"
    "wms_server_version"  "1.1.1"
    "wms_format"          "image/png"
    "wms_enable_request"  "*"
  END
  PROJECTION
     "init=epsg:900913"
  END
 
END

However if the service you are calling is secured it gets a whole load more complicated.  You might (as I did initially) think that all you have to do is change the map file to include the authorisation details such as

 
LAYER

NAME "MySecuredLayer"
DEBUG 5
  TYPE RASTER
  STATUS ON
  CONNECTION https://somedomain.com/WMSServer
  CONNECTIONTYPE WMS
  METADATA
    "wms_title"              "OSI"
    "wms_srs"             "EPSG:27700 EPSG:4326 EPSG:54004 EPSG:3978 EPSG:3857 EPSG:900913"
    "wms_name"            "1,2,3,4,5,6,7,8,9"
    "wms_server_version"  "1.1.1"
    "wms_format"          "image/jpeg"
    "wms_auth_username"   "myUserName"
    "wms_auth_password"   "myPassword"   
    "wms_auth_type"       "basic"
    "wms_enable_request"  "*"
  END
  PROJECTION
     "init=epsg:900913"
  END
 
END

Unfortunately that in itself doesn’t work (though the map file itself is correct).  You are very strongly recommended to make sure you have DEBUG setup – at least while testing.  So in addition to specifying the DEBUG level in the layer – you also need to include in the overall header of the MAP file a path to the log file

CONFIG "MS_ERRORFILE" "C:\\Mapserver\\ms_log.txt"
DEBUG 5

Running the above map file will create errors in the log file:

HTTP: request failed with curl error code 60 (SSL certificate problem, verify that the CA cert is OK. Details:
error:14090086:SSL routines:SSL3_GET_SERVER_CERTIFICATE:certificate verify failed) for https://somedomain.com/WMSServer…

(If you get an error:

 msHTTPExecuteRequests(): HTTP request error. HTTP: request failed with curl error code 1 (Protocol https not supported or disabled in libcurl)

Then you have a version of mapserver that is not compiled with the correct support and you need to change your version of mapserver (see below))

Basically the secure WMS call uses cURL libraries to do the call – and you have to register the certificate from the service you are calling.

What you need to do:

Make sure your version of mapserver supports https

Make sure you have a version of mapserver that supports https calls.  Tamas builds (http://www.gisinternals.com/sdk/) seem to include the required libraries. I used release-1500-gdal-1-10-0-mapserver-6-2-1.zip though I’m sure other versions would work as well

Get the security certificate

Use OpenSSL (whcih is included in the mapserver distribution in bin\curl)

openssl s_client –connect somedomain:443 > c:\Temp\cert.txt

If all goes well – the cert.txt file will look something like this – and you need the bits from -----BEGIN CERTIFICATE-----  to  -----END CERTIFICATE-----

CONNECTED(00000190)
---
Certificate chain
0 s:/OU=Domain Control Validated/OU=COMODO SSL Wildcard/CN=*.osi.ie
   i:/C=GB/ST=Greater Manchester/L=Salford/O=COMODO CA Limited/CN=COMODO SSL CA
1 s:/C=GB/ST=Greater Manchester/L=Salford/O=COMODO CA Limited/CN=COMODO SSL CA
   i:/C=SE/O=AddTrust AB/OU=AddTrust External TTP Network/CN=AddTrust External CA Root
---
Server certificate
-----BEGIN CERTIFICATE-----
MIIE5DCCA8ygAwIBAgIQOhI8AlhJBZcVfKqx4ruAnzANBgkqhkiG9w0BAQUFADBw
.
9WspBdI6j1t8l+ZZGe1PXBIBaHlpXuazvlC5FpU93bKt7zaEJ1iU5XjIwB/H4ys1
Hsch6g7YXnw=
-----END CERTIFICATE-----
subject=/OU=Domain Control Validated/OU=COMODO SSL Wildcard/CN=*.osi.ie
issuer=/C=GB/ST=Greater Manchester/L=Salford/O=COMODO CA Limited/CN=COMODO SSL CA
---
No client certificate CA names sent
---
SSL handshake has read 2675 bytes and written 536 bytes
---
New, TLSv1/SSLv3, Cipher is AES256-SHA
Server public key is 2048 bit
Secure Renegotiation IS NOT supported
Compression: NONE
Expansion: NONE
SSL-Session:
    Protocol  : TLSv1
    Cipher    : AES256-SHA
    Session-ID: 5F65EE6F4E78B564620101A5A706E7A259375B2277F430BC4578395780D14458
    Session-ID-ctx:
    Master-Key: 1C992E7CC5DFE1AF7B2D87295954D6A73A740E71BD81A2DB68546E61C055D0089FF1C91FB761F64CFE1BBD20B8411AE3
    Key-Arg   : None
    PSK identity: None
    PSK identity hint: None
    Start Time: 1391380021
    Timeout   : 300 (sec)
    Verify return code: 20 (unable to get local issuer certificate)
---

Add the security certificate to your local store

I copied the default file from bin\curl\curl-ca-bundle.crt into the main mapserver cgi-bin (or scripts) folder

Edit the curl-ca-bundle.crt file and add in your certificate (the bits from -----BEGIN CERTIFICATE-----  to  -----END CERTIFICATE-----)

Test your certificate

Check Curl.exe is in the same folder as your curl-ca-bundle.crt

1) I tested first with a call to getCapabilities:

curl -u MyUsername:MyPassword https://SomeDomain.com:443/WMSServer?REQUEST=GetCapabilities%26SERVICE=WMS

(Note I had to swap & for %26 0 i.e. url encode)

If this has worked you should see a normal getCapabilities text response

2) I then tested with an image call:

curl -u MyUsername:MyPassword https://SomeDomain.com:443/WMSServer?LAYERS=1%26REQUEST=GetMap%26SERVICE=WMS ……

If this works you should see a load of extended characters returned (i.e. an image displayed in the command window!)

Set a windows environment variable to point at the certificate store file

Previous tests have simply used the file in the same folder as the curl.exe, but when running for real mapserver must know where to find the certificate file.

Set the CURL_CA_BUNDLE environment variable to point at the certificate file

e.g. C:\Mapserver\cgi-bin\curl-ca-bundle.crt

Then you need to reboot the server to make sure the new environment setting is correctly used

Then test your map file

Now – at last – the map file which you thought should work originally should actually work and return you an image.

 
LAYER

NAME "MySecuredLayer"
DEBUG 5
  TYPE RASTER
  STATUS ON
  CONNECTION https://somedomain.com/WMSServer
  CONNECTIONTYPE WMS
  METADATA
    "wms_title"              "OSI"
    "wms_srs"             "EPSG:27700 EPSG:4326 EPSG:54004 EPSG:3978 EPSG:3857 EPSG:900913"
    "wms_name"            "1,2,3,4,5,6,7,8,9"
    "wms_server_version"  "1.1.1"
    "wms_format"          "image/jpeg"
    "wms_auth_username"   "myUserName"
    "wms_auth_password"   "myPassword"   
    "wms_auth_type"       "basic"
    "wms_enable_request"  "*"
  END
  PROJECTION
     "init=epsg:900913"
  END
 
END

If you look in the log file you should see:

Using CURL_CA_BUNDLE=C:\Mapserver\cgi-bin\curl-ca-bundle.crt

If not – check your environment variable setting

Other references

The above was put together from a variety of sources – including those below– useful additional refs include:

http://geographika.co.uk/setting-up-a-secure-cascading-wms-on-mapserver

http://www.mapserver.org/ogc/wxs_secure.html

http://trac.osgeo.org/mapserver/ticket/3070

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)

!table
!version 1000
!charset WindowsLatin1

Definition Table
  Type ODBC
begin_metadata
"\IsReadOnly" = "FALSE"
"\MapInfo" = ""
"\MapInfo\TableID" = "2e88f9e5-7d4c-4645-b824-4ef6b5102c9e"
"\DATALINK" = ""
"\DATALINK\Query" = "Select *  From ""myDatabase"".""dbo"".""myTableOrView"""
"\DATALINK\ConnectionString" = "Description=aDescription;DRIVER=SQL Server;SERVER=myServer;Trusted_Connection=Yes;APP=MapInfo Professional®;WSID=blah;DATABASE=myDatabase"
"\DATALINK\ToolKit" = "ODBC"
"\CACHE" = "OFF"
"\MBRSEARCH" = "ON"
end_metadata

lifted from here: http://stackoverflow.com/questions/2920948/selecting-a-good-sql-server-2008-spatial-index-with-large-polygons

Basically crate the SPROC listed at the bottom, then setup some example calls (3 or 4 different ones to minimise caching artefacts)

Then set it running and go on holiday …

Setup your parameters:

DECLARE @g1 VARCHAR(MAX)
SET @g1 = 'GEOMETRY::STGeomFromText(''POLYGON((252783.989267352 461095.507712082,260134.010732648 461095.507712082,260134.010732648 465068.492287918,252783.989267352 465068.492287918,252783.989267352 461095.507712082))'', 28992)'
DECLARE @g2 VARCHAR(MAX)
SET @g2 = 'GEOMETRY::STGeomFromText(''POLYGON((5580.146375 5340.100667,655613.186375 5340.100667,655613.186375 1219811.501845,5580.146375 1219811.501845,5580.146375 5340.100667))'',28992)'
exec sp_tune_spatial_index 'OBSVEG.OBS_Locations' ,'sidx_OBS_Locations_geom',2,8,@g1,@g2

Create Stored procedure:

N.B. Assumes geometry field is called geom
N.B. You MUST modify the SPROC to use extents for your SRID

USE [CMSI-OBSVEG]
GO
/****** Object:  StoredProcedure [dbo].[sp_tune_spatial_index]    Script Date: 10/01/2013 11:35:27 ******/
SET ANSI_NULLS ON
GO
SET QUOTED_IDENTIFIER ON
GO

ALTER PROCEDURE [dbo].[sp_tune_spatial_index]
(
  @tabnm                VARCHAR(MAX),    -- This parameter stores the name of the spatial table for which you are tuning the index
    @idxnm                VARCHAR(MAX),    -- This parameter stores the name of the spatial index of the named table
    @min_cells_per_obj    INT,            -- Minimum Cells Per Object to test on. Suggested to start at 2.
    @max_cells_per_obj    INT,            -- Maximum Cells Per Object to test on.
    
    /*    The test requires two geometry instances to use in test query 1 and 2.
        The first one should cover the area of default extent. The second should
        cover an area roughly the size of the area shown when zoomed in, panning
        around. It is required that the variable store a string that will create
        the geometry instance since this will be done within the procedure and 
        cannot be a variable of type: GEOMETRY. The SRID of these instances must
        match that of the table you are testing. */
    @testgeom1            VARCHAR(MAX),    -- This parameter stores the first geometry instance creation string that will be used in the test
    @testgeom2            VARCHAR(MAX)    -- This parameter stores the second geometry instance creation string that will be used in the test
    
)

AS

SET NOCOUNT ON;

/*    Prior to running this procedure, two tables are required. These tables are 
    created here to prepare for running the procedure.    */

PRINT 'Checking for required tables...'
IF EXISTS(SELECT 1 FROM sysobjects WHERE name IN ('cell_opt_perm', 'spat_idx_test_result'))
    BEGIN
        PRINT '... The "cell_opt_perm" and "spat_idx_test_result" tables exist.'
    END
ELSE
BEGIN
    PRINT '... Creating "cell_opt_perm" and "spat_idx_test_result" tables.'
    CREATE TABLE cell_opt_perm(
        [perm_id] [smallint] NOT NULL,
        [permutation] [nvarchar](4) NOT NULL,
        [level1] [nvarchar](6) NOT NULL,
        [level2] [nvarchar](6) NOT NULL,
        [level3] [nvarchar](6) NOT NULL,
        [level4] [nvarchar](6) NOT NULL
    )

    INSERT INTO cell_opt_perm ([perm_id], [permutation], [level1], [level2], [level3], [level4])
    VALUES (1,'LLLL','LOW','LOW','LOW','LOW'),
        (2,'LLLM','LOW','LOW','LOW','MEDIUM'),
        (3,'LLLH','LOW','LOW','LOW','HIGH'),
        (4,'LLML','LOW','LOW','MEDIUM','LOW'),
        (5,'LLMM','LOW','LOW','MEDIUM','MEDIUM'),
        (6,'LLMH','LOW','LOW','MEDIUM','HIGH'),
        (7,'LLHL','LOW','LOW','HIGH','LOW'),
        (8,'LLHM','LOW','LOW','HIGH','MEDIUM'),
        (9,'LLHH','LOW','LOW','HIGH','HIGH'),
        (10,'LMLL','LOW','MEDIUM','LOW','LOW'),
        (11,'LMLM','LOW','MEDIUM','LOW','MEDIUM'),
        (12,'LMLH','LOW','MEDIUM','LOW','HIGH'),
        (13,'LMML','LOW','MEDIUM','MEDIUM','LOW'),
        (14,'LMMM','LOW','MEDIUM','MEDIUM','MEDIUM'),
        (15,'LMMH','LOW','MEDIUM','MEDIUM','HIGH'),
        (16,'LMHL','LOW','MEDIUM','HIGH','LOW'),
        (17,'LMHM','LOW','MEDIUM','HIGH','MEDIUM'),
        (18,'LMHH','LOW','MEDIUM','HIGH','HIGH'),
        (19,'LHLL','LOW','HIGH','LOW','LOW'),
        (20,'LHLM','LOW','HIGH','LOW','MEDIUM'),
        (21,'LHLH','LOW','HIGH','LOW','HIGH'),
        (22,'LHML','LOW','HIGH','MEDIUM','LOW'),
        (23,'LHMM','LOW','HIGH','MEDIUM','MEDIUM'),
        (24,'LHMH','LOW','HIGH','MEDIUM','HIGH'),
        (25,'LHHL','LOW','HIGH','HIGH','LOW'),
        (26,'LHHM','LOW','HIGH','HIGH','MEDIUM'),
        (27,'LHHH','LOW','HIGH','HIGH','HIGH'),
        (28,'MLLL','MEDIUM','LOW','LOW','LOW'),
        (29,'MLLM','MEDIUM','LOW','LOW','MEDIUM'),
        (30,'MLLH','MEDIUM','LOW','LOW','HIGH'),
        (31,'MLML','MEDIUM','LOW','MEDIUM','LOW'),
        (32,'MLMM','MEDIUM','LOW','MEDIUM','MEDIUM'),
        (33,'MLMH','MEDIUM','LOW','MEDIUM','HIGH'),
        (34,'MLHL','MEDIUM','LOW','HIGH','LOW'),
        (35,'MLHM','MEDIUM','LOW','HIGH','MEDIUM'),
        (36,'MLHH','MEDIUM','LOW','HIGH','HIGH'),
        (37,'MMLL','MEDIUM','MEDIUM','LOW','LOW'),
        (38,'MMLM','MEDIUM','MEDIUM','LOW','MEDIUM'),
        (39,'MMLH','MEDIUM','MEDIUM','LOW','HIGH'),
        (40,'MMML','MEDIUM','MEDIUM','MEDIUM','LOW'),
        (41,'MMMM','MEDIUM','MEDIUM','MEDIUM','MEDIUM'),
        (42,'MMMH','MEDIUM','MEDIUM','MEDIUM','HIGH'),
        (43,'MMHL','MEDIUM','MEDIUM','HIGH','LOW'),
        (44,'MMHM','MEDIUM','MEDIUM','HIGH','MEDIUM'),
        (45,'MMHH','MEDIUM','MEDIUM','HIGH','HIGH'),
        (46,'MHLL','MEDIUM','HIGH','LOW','LOW'),
        (47,'MHLM','MEDIUM','HIGH','LOW','MEDIUM'),
        (48,'MHLH','MEDIUM','HIGH','LOW','HIGH'),
        (49,'MHML','MEDIUM','HIGH','MEDIUM','LOW'),
        (50,'MHMM','MEDIUM','HIGH','MEDIUM','MEDIUM'),
        (51,'MHMH','MEDIUM','HIGH','MEDIUM','HIGH'),
        (52,'MHHL','MEDIUM','HIGH','HIGH','LOW'),
        (53,'MHHM','MEDIUM','HIGH','HIGH','MEDIUM'),
        (54,'MHHH','MEDIUM','HIGH','HIGH','HIGH'),
        (55,'HLLL','HIGH','LOW','LOW','LOW'),
        (56,'HLLM','HIGH','LOW','LOW','MEDIUM'),
        (57,'HLLH','HIGH','LOW','LOW','HIGH'),
        (58,'HLML','HIGH','LOW','MEDIUM','LOW'),
        (59,'HLMM','HIGH','LOW','MEDIUM','MEDIUM'),
        (60,'HLMH','HIGH','LOW','MEDIUM','HIGH'),
        (61,'HLHL','HIGH','LOW','HIGH','LOW'),
        (62,'HLHM','HIGH','LOW','HIGH','MEDIUM'),
        (63,'HLHH','HIGH','LOW','HIGH','HIGH'),
        (64,'HMLL','HIGH','MEDIUM','LOW','LOW'),
        (65,'HMLM','HIGH','MEDIUM','LOW','MEDIUM'),
        (66,'HMLH','HIGH','MEDIUM','LOW','HIGH'),
        (67,'HMML','HIGH','MEDIUM','MEDIUM','LOW'),
        (68,'HMMM','HIGH','MEDIUM','MEDIUM','MEDIUM'),
        (69,'HMMH','HIGH','MEDIUM','MEDIUM','HIGH'),
        (70,'HMHL','HIGH','MEDIUM','HIGH','LOW'),
        (71,'HMHM','HIGH','MEDIUM','HIGH','MEDIUM'),
        (72,'HMHH','HIGH','MEDIUM','HIGH','HIGH'),
        (73,'HHLL','HIGH','HIGH','LOW','LOW'),
        (74,'HHLM','HIGH','HIGH','LOW','MEDIUM'),
        (75,'HHLH','HIGH','HIGH','LOW','HIGH'),
        (76,'HHML','HIGH','HIGH','MEDIUM','LOW'),
        (77,'HHMM','HIGH','HIGH','MEDIUM','MEDIUM'),
        (78,'HHMH','HIGH','HIGH','MEDIUM','HIGH'),
        (79,'HHHL','HIGH','HIGH','HIGH','LOW'),
        (80,'HHHM','HIGH','HIGH','HIGH','MEDIUM'),
        (81,'HHHH','HIGH','HIGH','HIGH','HIGH')
    
    CREATE TABLE spat_idx_test_result(
        [perm_id] [int] NOT NULL,
        [num_cells] [int] NOT NULL,
        [permut] [nvarchar](4) NOT NULL,
        [g1t1] [bigint] NULL,
        [g1t2] [bigint] NULL,
        [g1t3] [bigint] NULL,
        [g1t4] [bigint] NULL,
        [g2t1] [bigint] NULL,
        [g2t2] [bigint] NULL,
        [g2t3] [bigint] NULL,
        [g2t4] [bigint] NULL
    )
    
    INSERT INTO dbo.spat_idx_test_result
    VALUES (0,16,0,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)
END


/*    delete all rows from "spat_idx_test_result" table. This makes it ready to stuff in new results.
    !!!WARNING!!! if your test was interupted, the table will be cleared out and the test will
    begin from the beginning. You could try to modify this to start where you left off but
    I didn't have time and this worked well enough for me. */
DELETE FROM spat_idx_test_result
    WHERE perm_id != 0

/* set up counters */
DECLARE @a1 INT
DECLARE @a2 INT
DECLARE @a3 INT
DECLARE @a4 INT

/* set up variables to hold high/medium/low values and permutation to use in rebuilding
   the spatial index and recording stats */
DECLARE @lev1 VARCHAR(6)
DECLARE @lev2 VARCHAR(6)
DECLARE @lev3 VARCHAR(6)
DECLARE @lev4 VARCHAR(6)
DECLARE @permut VARCHAR(6)
DECLARE @num_cell VARCHAR(4)
DECLARE @time_str VARCHAR(20)
DECLARE @perm_id VARCHAR(20)

/* create variables to hold timestamps for beginning and ending of test queries */
DECLARE @start_t DATETIME
DECLARE @end_t DATETIME
DECLARE @elapse_t INT

/* begin looping through cell option permutations */
SET @a1 = @min_cells_per_obj
WHILE @a1 <= @max_cells_per_obj
    BEGIN
        SET @a2 = 1
        PRINT 'Started Testing for ' +CAST(@a1 AS VARCHAR(10)) +' cells per object'
        WHILE @a2 < 82
            BEGIN
                SELECT @lev1 = level1, @lev2 = level2, @lev3 = level3, @lev4 = level4 FROM cell_opt_perm WHERE perm_id = @a2
                SET @permut = '''' +(SELECT permutation FROM cell_opt_perm WHERE perm_id = @a2) +''''
                EXEC
                    ('
                        CREATE SPATIAL INDEX ' +@idxnm +' ON ' +@tabnm +' 
                        (
                            [geom]
                        )
                        USING  GEOMETRY_GRID 
                        WITH
                            (
                                BOUNDING_BOX =(-100, -100, 300000, 650000),
                                GRIDS =(LEVEL_1 = ' +@lev1 +' ,LEVEL_2 = ' +@lev2 +' ,LEVEL_3 = ' +@lev3 +' ,LEVEL_4 = ' +@lev4 +' ), 
                                CELLS_PER_OBJECT = ' +@a1 +' ,
                                PAD_INDEX  = OFF,
                                SORT_IN_TEMPDB = OFF,
                                DROP_EXISTING = ON,
                                ALLOW_ROW_LOCKS  = ON,
                                ALLOW_PAGE_LOCKS  = ON,
                                FILLFACTOR = 100
                            )
                        ON [PRIMARY]'
                    )
                PRINT 'Re-built index to ' +@permut
                SET @a3 = 1
                SET @a4 = 1
                WHILE @a3 < 5
                    BEGIN
                        SET @start_t = GETDATE()
                        EXEC
                            (
                                'CREATE TABLE #tmp_tab (shp GEOMETRY)
                                DECLARE @g1 GEOMETRY
                                SET @g1 = ' +@testgeom1 +'
                                INSERT #tmp_tab (shp)
                                    SELECT
                                        r.geom AS shp
                                    FROM
                                        ' +@tabnm +' r
                                    WHERE
                                        r.geom.STIntersects(@g1) = 1
                                DROP TABLE #tmp_tab'
                            )
                        SET @end_t = GETDATE()
                        SET @elapse_t = (SELECT DATEDIFF(MS, @start_t, @end_t))
                        SET @num_cell = CAST(@a1 AS VARCHAR(6))
                        SET @time_str = CAST(@elapse_t AS VARCHAR(20))
                        IF @a3 = 1
                            BEGIN
                                IF (SELECT TOP 1 perm_id FROM spat_idx_test_result) IS NULL
                                    BEGIN
                                        SET @perm_id = 1
                                    END
                                ELSE
                                    BEGIN
                                        SET @perm_id = CAST((SELECT MAX(perm_id+1) FROM spat_idx_test_result) AS VARCHAR(20))
                                    END
                                EXEC
                                    (
                                        'INSERT INTO spat_idx_test_result (perm_id, num_cells, permut, g1t' +@a3 +')
                                        VALUES (' +@perm_id +', ' +@num_cell +', ' +@permut +', ' +@time_str +')'
                                    )
                            END
                        ELSE
                            EXEC
                                (
                                    'UPDATE spat_idx_test_result
                                    SET
                                    num_cells = ' +@num_cell +',
                                    permut = ' +@permut +',
                                    g1t' +@a3 +' = ' +@time_str +'
                                    WHERE perm_id = ' +@perm_id
                                )
                        SET @a3 = @a3 + 1
                    END
                WHILE @a4 < 5
                    BEGIN
                        SET @start_t = GETDATE()
                        EXEC
                            (
                                'CREATE TABLE #tmp_tab (shp GEOMETRY) 
                                DECLARE @g2 GEOMETRY
                                SET @g2 = ' +@testgeom2 +'
                                INSERT #tmp_tab (shp)
                                    SELECT
                                        r.geom AS shp
                                    FROM
                                        ' +@tabnm +' r
                                    WHERE
                                        r.geom.STIntersects(@g2) = 1
                                DROP TABLE #tmp_tab'
                            )
                        SET @end_t = GETDATE()
                        SET @elapse_t = (SELECT DATEDIFF(MS, @start_t, @end_t))
                        SET @num_cell = CAST(@a1 AS VARCHAR(6))
                        SET @time_str = CAST(@elapse_t AS VARCHAR(20))
                        EXEC
                            (
                                'UPDATE spat_idx_test_result
                                SET
                                num_cells = ' +@num_cell +',
                                permut = ' +@permut +',
                                g2t' +@a4 +' = ' +@time_str +'
                                WHERE perm_id = ' +@perm_id
                            )
                        SET @a4 = @a4 + 1
                    END
                SET @a2 = @a2 + 1
            END
        SET @a1 = @a1 + 1
    END
PRINT 'Testing of ' +@tabnm +' spatial index: ' +@idxnm +' is complete!'

1. the WMS request should include this on the KVP URL: &tiled=true
2. the client should construct its tile requests so that they match the server tiles – this is not the place for elaboration, but essentially this means getting the origin and tile size correct
3. the GeoServer must have “Enable direct integration with GeoServer WMS” switched on.

The setting in #3 is found under Tile Caching > Caching Defaults in the main GeoServer menu. Set it once, and forget it. Or so I thought.

How do we know if the WMS is using its cache? Insect the http requests in Fiddler (or similar) and check the headers of the responses.

If the WMS is using its cache you will see lines in the headers like this

geowebcache-cache-result: HIT

along with various other geowebcache-… values

If it knows it has a cache, but something is wrong with the WMS requests (item #2 above), you might see headers like this:

geowebcache-miss-reason: request does not align to grid(s) 'OSGB27700'
geowebcache-cache-result: MISS

And if the WMS is unaware it even has a cache of tiles at all, you will see nothing about geowebcache at all in the headers. And of course, poor WMS performance may give you a clue something is not working as you hoped.

Well I had set everything up lovely and my cache was working well, giving super fast performance of my Ordnance Survey mapping once the initial requests had populated the cache for an area.

But the very next day… performance was rubbish, and on inspecting the http headers, it turned out the WMS had forgotten about its cache.

In the GeoServer admin interface, the “Enable direct integration with GeoServer WMS” was now set to false. Grrr – I wondered how I could stupidly have switched that off again. Ah well, switch it on again and move on.

A few hours later… same again. This time I started suspecting colleagues, but there was no evidence, or motive, yet…

A couple more times, and I realised it was reverting after any re-start of GeoServer. I had been restarting every now and then, to make sure changed wrapper log settings etc. stuck, but it turned out this had been un-sticking my essential cache setting.

The solution turned out to be to edit the gwc-gs.xml file, found in the GeoServer data directory, and set this:

  <directWMSIntegrationEnabled>true</directWMSIntegrationEnabled>

In my data dir, there was a gwc-gs.xml.tmp file, and I can only assume GeoServer had always failed to update the main XML file for some reason, used the .tmp version for the duration of the session, then reverted to the main one upon a re-start.

Now my cached WMS layers run nice and fast, and GeoServer doesn’t forget the setting.

In order to resolve this, MapInfo introduced a 1 mm accuracy projection, which can be found under Projection... when you create a table or save a copy of an existing table. However, this is not a default option, so even when saving a copy of data already in the 1 mm precision projection MapInfo will automatically default to the standard projection. This makes it very easy to accidentally reduce the quality of your data if you forget to actively select the 1 mm precision projection.

Fortunately there is a relatively simple way to ensure that MapInfo always defaults to the 1 mm projection:

• Navigate to the MapInfo folder under program files (e.g. C:\Program Files (x86)\MapInfo\Professional).
• Locate the MAPINFOW.PRJ file.
• Backup the MAPINFOW.PRJ file, in case you make a mistake.
• Open the MAPINFOW.PRJ file in a text editor (e.g. notepad).
• Search for ‘British’, which will find a section of text that looks like this:

"--- British Coordinate Systems ---"
"British National Grid\p27700", 8, 79, 7, -2, 49, 0.9996012717, 400000, -100000
"British National Grid (1 mm accuracy)", 2008, 79, 7, -2, 49,0.9996012717, 400000, -100000, 0, 0, 2000000, 2000000
"British OSGRS80 Grid", 8, 33, 7, -2, 49, 0.9996012717, 400000, -100000

• Move the third line of this text up one so that it becomes the second:

"--- British Coordinate Systems ---"
"British National Grid (1 mm accuracy)", 2008, 79, 7, -2, 49,0.9996012717, 400000, -100000, 0, 0, 2000000, 2000000
"British National Grid\p27700", 8, 79, 7, -2, 49, 0.9996012717, 400000, -100000
"British OSGRS80 Grid", 8, 33, 7, -2, 49, 0.9996012717, 400000, -100000

This change moves the 1 mm accuracy projection to the top of the Category Members selector under the British Coordinate Systems Category in the Choose Projection dialogue box of MapInfo, automatically making it the default option (as shown in the picture). If you want to use the standard projection then you just need to select 'British National Grid [EPSG: 27700]' from the dialogue instead.

This is a useful script, using CTEs (common table expressions) to calculate the extent of the geometry in a SQL Server table. Replace “geom” with the name of your geometry column and “MyTable” with the name of your table.

Note that this script depends on STEnvelope returning a rectangular extent with its first point at the SW corner and its third in the NE corner. Tests from the author of the blog above confirm the SQL Server currently does this.

with
  cte_envelope as
  (
    select geom.STEnvelope() as envelope from MyTable
  ),
  cte_corner as
  (
    select envelope.STPointN(1) as point from cte_envelope
    union all
    select envelope.STPointN(3) from cte_envelope
  )
select min(point.STX) as min_x, min(point.STY) as min_y, max(point.STX) as max_x, max(point.STY) as max_y
from cte_corner

Output quality and performance

	GeoServer	MapServer
Response time	0.8 to 1.0 seconds with USE_JAI_IMAGEREAD=true 0.6 to 1.3 seconds with USE_JAI_IMAGEREAD=false	0.4 to 0.6 seconds
Image size	63,574 bytes	78,327 bytes
The map
Quality	Adequate, but quite fractured	Lovely
Comments	A significantly slower response than MapServer, and a worse image. I couldn’t find any difference between the different re-sampling methods, perhaps suggesting that the settings were not taking effect (I’ve seen GeoServer GUI issues like this a couple of times now, requiring diving into the XML config files). Later – yes the settings were present in the XML config, but seem to make no difference to the output.	Faster and much nicer image. If we turn off bilinear re-sampling, then the response time improves to around 0.3 seconds and the image deteriorates slightly as shown below (becoming similar in quality to the GeoServer map).

Oh dear this is not looking good for GeoServer. But is it the whole story? Next I compared them in a more realistic scenario, as base maps in OpenLayers, trying both single tile and multi-tile modes. I configured the map to view the 1:50,000 mapping at a scale where no re-sampling could occur.

	GeoServer	MapServer
Response time multi-tile (30 tiles each 256x256 pixels)	0.2 to 1.2 seconds, but on average a shade faster	0.3 to 0.7 seconds
Response time single-tile (1725 x 1173 pixels)	4 to 6 seconds	3 to 5 seconds
Image size (typical 256x256 tile)	30,447 bytes	22,674 seconds
Image size (1725 x 1173 pixels)	646,681 bytes	690,494 bytes
The map (256x256 tile)
Quality	Very good	Very good (pretty much identical)

OK so now it is looking better for GeoServer; when there is no re-sampling it can return identical images to MapServer, though slightly larger, and performance is pretty similar and possibly quicker.

Output quality and performance when re-projecting

Now time to see how they cope with serving the maps into a different coordinate system, and in particular the global spherical Mercator used in Google and Bing Maps.

	GeoServer	MapServer
Response time multi-tile (30 tiles each 256x256 pixels)	0.2 to 1.2 seconds	0.1 to 0.5 seconds
Response time single-tile (1725 x 1173 pixels)	10 seconds	6 seconds
Image size (typical 256x256 tile)	29,169 bytes	37,261 seconds
Image size (1725 x 1173 pixels)	1,041,098 bytes	982,050 bytes
The map (256x256 tile)
Quality	Pretty horrible	Very good

With default settings, GeoServer produces horrible output when re-projecting the raster maps. But this does need more investigation when I found out how to make it use different re-sampling algorithms.

Performance under load

Now time to see how they perform under load using multi-mechanize. I chose to compare MapServer, GeoServer and a tilecache at returning a 256x256 tile in EPSG:900913 (i.e. re-projecting on-the-fly). Each test ran for 100 seconds, starting with 10 virtual users, increasing by 10 every 10 seconds. The tilecache image was of course cached and simply being retrieved from disk. All tests were using the same web server, with files on the same drive.

Note: the vertical scale differs in each graph to fit the observed values.

	GeoServer	MapServer	Tilecache
Transactions	381	556	2188
Errors	57	0	0
Average response time	6.1	6.5	2.3
95% response time	26.8	17.6	5.1
Response graph
Transactions per second
Comments	Completely broken after about 80 seconds – had to re-start GeoServer	No errors, but significantly disrupted performance	Slight slow-down under high load, but very reassuring performance

GeoServer clearly could not stand up to more than about 30 virtual users when re-projecting the maps, and ultimately the service stalled completely requiring a re-start.

I tried the tests again requesting the maps in OSGB so that the service did not involve re-projection. I also reduced the number of virtual users to 50. The results were quite surprising…

	GeoServer	MapServer	Tilecache
Transactions	2018	1114	2197
Errors	0	0	0
Average response time	1.3	2.4	1.1
95% response time	2.4	9.1	2.5
Response graph
Transactions per second
Comments	Impeccable!	Much better than when re-projecting, but only half as good as GeoServer.	Only fractionally better than GeoServer

This is a hugely impressive result for GeoServer. My guess is that it is employing server-side in-memory caching, and because all my requests were identical the responses were very fast. MapServer on the other hand has to start a CGI process of each request and presumably cannot benefit from an in-memory cache.

When re-projecting on the other hand, perhaps it does not use a cache. And because all the GeoServer operations are running as one process, it is restricted to one virtual processor on our server whereas each MapServer exe can grab one for itself. Certainly when running these tests the server was running at about 40-80% CPU with MapServer, compared to 7% with GeoServer.

I ran this test again with the number of virtual users doubled back to 100, and GeoServer came through almost unscathed, this time processing 2072 transactions though with 3 errors.

	GeoServer	MapServer	Tilecache
Transactions	2072	907	2234
Errors	3	0	0
Average response time	2.4	5.9	2.3
95% response time	5.9	25.6	4.8
Response graph
Transactions per second
Comments	Wow	Suffering, but no errors	The best as expected

I suspect this test may not predict real-world behaviour with lots of users requesting different WMS maps, where I suspect GeoServer would lose its advantage, but it is interesting anyway.

NB I had GeoServer running with AllowMultithreading=true

Summary

If the maps need to be re-projected then use MapServer – both performance and output quality are vastly superior. In a site that is going to get any significant load it is unwise to re-project raster mapping on-the-fly anyway, and it should definitely be cached. In this case the WMS performance doesn’t matter so much but output quality is paramount, so MapServer it is.  If I find a way of making GeoServer produce nice looking output I’ll come back and revise this post!

If no re-projection is required, then simple observation has the two roughly equal under light load. The maps looked almost identical, and performance was close. GeoServer had the advantage for me in that it didn’t suffer from MapServer’s tendency to choke on some requests; in practice this meant I went over my entire OSGB cache again with a GeoServer WMS as the source, to plug the gaps left by the MapServer service.

Under heavy (albeit contrived) load, the crucial factor was whether or not there was any re-projection involved. Without re-projection GeoServer produced astonishing performance, thrashing MapServer and nearly matching a tilecache. Introduce re-projection however, and GeoServer collapses to the extent that it can crash entirely – not good when it is running all your mapping services. MapServer never broke, and despite some gaps and slow responses, it managed to struggle through all tests without a single error.

So out of this experience I think we need to be using both, depending on the purpose.

Next to test performance on vector layers, including file based GIS formats and spatial databases…

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.

Laine Skinner from ESDM presented a paper about a mobile application we have written for our HBSMR product, which is used by Historic Environment Records to manage heritage inventories and to support casework.

The presentation described how we created a mobile application over a short-targeted project in January 2011. The project itself was completed in a RAD/agile manner and highlights a successful way of working with a technology which itself is also undergoing huge change every few months.

The presentation is available online here.

Unfortunately there is no online content from the conference, although speaker details are available from the AGI and speakers are usually more than happy to discuss their work and forward details onward.

Laine’s presentation was titled:

GIS Cloud: Money for old rope or a new direction?

Its available as a prezi (cloud based presentation service) here. The presentation provides a short description about what the ‘GIS Cloud’ is and why it might be useful to a whole range of organizations. We discuss some practical uses of a real cloud implementation on Amazons AWS product.

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…

• If you have our ArcGIS Maplink product (or ArcSquirrel) then loading data in using this into a temporary SQL table is probably best as it will do a reasonable amount of checking.
• MapInfo has its own SQL loader tool ("Easy Loader" - find it in the Tool Manager)
• If you are just using open source you can use OGR2OGR

ogr2ogr -overwrite -f "MSSQLSpatial" "MSSQL:server=MyServer;database=MyDatabase;Trusted_Connection=yes;" "C:\Temp\MyShapeFile.shp" -a_srs "EPSG:27700"

Once the data is in SQL there can still be issues:

Incorrect SRID

• Make sure you specified the correct SRID when the data was loaded (if you didn’t specify it – it may well be 0). You can check this with the SQL: SELECT geom.STSrid from SitesGeom Where geom.STSrid <> 27700

Invalid ring

• You may also have invalid ring geometries – a good post on this is http://blogs.msdn.com/b/edkatibah/archive/2008/08/19/working-with-invalid-data-and-the-sql-server-2008-geography-data-type-part-1b.aspx but the key SQL command is:

UPDATE foo SET geom = geom.STUnion(geom.STStartPoint());

Invalid geometries

UPDATE foo SET geom = geom.MakeValid()

SRID not imported

Import geometry data into temp SQL table (e.g. foo) – with tempGeom column

Then try to populate a geometry column from the WKB of the validated column:

ALTER TABLE foo ADD geomFinal geometry;

SET geomFinal = geometry::STGeomFromWKB(geomTemp.STAsBinary(), 27700)

http://www.ordnancesurvey.co.uk/oswebsite/products/geo-referencing.html

These are available as *.tfw "world files" for ArcGIS^®, or *.tab files for MapInfo Professional^®

We have written a custom MapBasic app for this purpose: for ESDM staff this is located in 3rdPartyApps\TABs_from_ECW [this being a pretty dumb place for it, as it is not a 3rdPartyApp, Ed]. Instructions are in the text file in the zip.

If you are not Exegesis staff, and find yourself in this situation, please contact us for help.

I think it uses ogr2ogr under the bonnet.

http://mygeodata.eu/apps/converter/index_EN.html