PostGIS Spatial Database Engine UMN Mapserver Boston Geographic Information Systems    PostGreSQL Object Relational Database Management System
Home   Consulting Services  About Boston GIS   Boston GIS Blog  Postgres OnLine Journal

Table Of Contents

UMN Mapserver Tutorials

How to Use different kinds of datasources in UMN Mapserver layers
Using Mapserver as a WMS Client.

PostGIS and PostgreSQL Tutorials

Solving the Nearest Neighbor Problem in PostGIS
PostGIS Nearest Neighbor: A Generic Solution - Much Faster than Previous Solution
Part 1 - PostGIS and SharpMap in ASP.NET 2.0 using VB.NET: Compiling SharpMap with PostGIS
Part 2 - PostGIS and SharpMap in ASP.NET 2.0 using VB.NET: Displaying the Maps
Part 1: Getting Started With PostGIS: An almost Idiot's Guide
Part 2: Introduction to Spatial Queries and SFSQL
Part 3: PostGIS Loading Data from Non-Spatial Sources
PLR Part 1: Up and Running with PL/R (PLR) in PostgreSQL: An almost Idiot's Guide
PLR Part 2: PL/R and PostGIS
PLR Part 3: PL/R and Geospatial Data Abstraction Library (GDAL) RGDAL

Miscellaneous Tutorials/Cheatsheets/Examples

OGR2OGR Cheatsheet
Using OpenLayers: Part 1
Using OpenLayers: Part 2

PostGIS Code Snippets

PostGIS ST_Dump, Dump
Extent Expand Buffer Distance: PostGIS - ST_Extent, Expand, ST_Buffer, ST_Distance
PostGIS GeomFromText
Intersects Intersection: PostGIS - ST_Intersects, ST_Intersection
PostGIS MakeLine ST_MakeLine Examples
PostGIS MakePoint
PointFromText, LineFromText, ST_PointFromText, ST_LineFromText OGC functions - PostGIS
PostGIS Simplify
Summary (pre 1.3.1 name), ST_Summary (+1.3.1)
PostGIS: ST_Translate

UMN Mapserver Tutorials

How to Use different kinds of datasources in UMN Mapserver layers

One of the very great things about the UMN Mapserver Web system is that it can support numerous kinds of datasources. In this brief excerpt we will provide examples of how to specify the more common data sources used for layers. The examples below are for Mapserver 4.6, but for the most part are applicable to lower versions.

File locations for file based datasources such as ESRI Shape and MapInfo tab files are defined relative to the SHAPEPATH attribute or as absolute paths. For example the beginning declaration of your .map file might look something like the below


MAP 
    #
    # Start of map file
    #

    NAME MYMAP
    
    EXTENT  732193.725550 2904132.702662 799614.090681 2971466.288170 
    
    SIZE 500 500
    SHAPEPATH "c:\mydata\"
    :
    :

ESRI Shapefile

The most common kind of data used in UMN Mapserver is the ESRI shapefile which has a .shp extension. For this kind of datasource you simply specify the location of the file without even specifying the extension. Below is a sample declaration of a polygon layer that uses a shape file

    LAYER
        NAME buildings
        TYPE POLYGON
        STATUS DEFAULT
        DATA buildings
        PROJECTION
            "init=epsg:2249"
        END
        CLASS
            OUTLINECOLOR 10 10 10
        END
    END

MapInfo Tab Files

Many datasources are available to mapserver via the GDAL OGR driver. Map Info is one of those datasources. Below example is what a mapinfo layer definition looks like. Note the tab file specified should be placed in the folder denoted by SHAPEPATH at top of map file



    LAYER
        NAME buildings
        STATUS DEFAULT
        MINSCALE 7000
        CONNECTIONTYPE OGR
        CONNECTION "buildings.tab" 
        TYPE POLYGON
        PROJECTION
            "init=epsg:2249"  
        END
        # -- MapInfo has projection information built in the tab file 
        # -- so you can often auto read this information with the below
        #PROJECTION
        # AUTO
        #END
        CLASS
            OUTLINECOLOR 10 10 10
        END
    END

PostGIS Layer

Mapserver has a custom driver for the PostGIS spatial database. In order to use this, your mapserver cgi or mapscript must be compiled with the PostGIS driver. Below is what a postgis mapserver layer looks like.



    LAYER
      CONNECTIONTYPE postgis
      NAME "buildings"
      CONNECTION "user=dbuser dbname=mydb host=myserver"
      # the_geom column is the name of a spatial geometry field in the table buildings
      DATA "the_geom from buildings"
      STATUS DEFAULT
      TYPE POLYGON
      # Note if you use a filter statement - this is basically like a where clause of the sql statement
      FILTER "storyhg > 2"
      CLASS
            OUTLINECOLOR 10 10 10
      END
    END

More complex PostGIS layer

LAYER
    NAME "projects"
    CONNECTIONTYPE postgis
    CONNECTION "user=myloginuser dbname=mydbname host=mydbhost password=mypass" 
    DATA "the_geom FROM (SELECT a.projid, a.projname, a.projtype, a.projyear, a.pid, parc.the_geom 
			FROM projects a INNER JOIN parcels parc ON a.parcel_id = parc.pid 
                         WHERE a.projyear = 2007) as foo USING UNIQUE projid USING SRID=2249" 
    STATUS OFF
    TYPE POLYGON
    CLASS
        NAME "Business Projects"
        EXPRESSION ('[projtype]' = 'Business')
        STYLE
            OUTLINECOLOR 204 153 51
            WIDTH 3
        END
    END
    CLASS
        NAME "Community Projects"
        EXPRESSION ('[projtype]' = 'Community')
        STYLE
            OUTLINECOLOR 204 0 0
            WIDTH 3
        END
    END

    PROJECTION
        "init=epsg:2249" 
    END
    METADATA
     "wms_title" "Projects"
     "wfs_title" "Projects"
      gml_include_items "all"
      wms_include_items "all"
    END
    DUMP TRUE
    TOLERANCE 10
END

WMS Layer

Mapserver has the ability to act as a WMS Server as well as a WMS Client. The WMS Client capabilities are accessed by defining WMS layers that connect to WMS servers. Below is an example of a WMS layer using the Microsoft Terraservices WMS Server.



    LAYER
          NAME "msterraservicedoq"
          TYPE RASTER
          STATUS DEFAULT
          CONNECTION "http://terraservice.net/ogcmap.ashx?"
          CONNECTIONTYPE WMS
          MINSCALE 3000
          MAXSCALE 20000
          #DEBUG ON
          METADATA
            "wms_srs"             "EPSG:26919"
            "wms_name"            "doq"
            "wms_server_version"  "1.1.1"
            "wms_format"          "image/jpeg"
            "wms_style"         "UTMGrid_Cyan"
            "wms_latlonboundingbox" "-71.19 42.23 -71 42.40"
          END
     END


Using Mapserver as a WMS Client.

Example mapserver map that calls microsoft terraservice WMS.

In this example we show how to use Mapserver as a WMS client by utilizing Microsoft's Terra Service WMS server. For more details about Microsft's OGC WMS check out the GetCapabilities of Microsoft Terraservice.

download

PostGIS and PostgreSQL Tutorials

Solving the Nearest Neighbor Problem in PostGIS

A common problem encountered in GIS is the Nearest Neighbor problem. In a nutshell the problem is to find the x number of nearest neighbors given a geometry and n geometries of data.

The nearest neighbor crops up in other disciplines as well, except in other disciplines the units of measurement are not spatial distance but instead some sort of matrix distance. For this particular talk, I will focus on the GIS nature of the problem and specifically solving it in PostGIS.

The easiest and most naive way to do this is to do an exhaustive search of the whole space comparing the distances against the geometry of interest and sorting them. This solution while intuitive and easy to write is the slowest way of doing this.

If you were to write such a simple solution in PostgreSQL, then it would take something of the form shown below. The below gives you the nearest 5 neighbors from a reference row with gid = 1.

SELECT g1.gid As gref_gid, g1.description As gref_description, 
g2.gid As gnn_gid, g2.description As gnn_description 
    FROM sometable As g1, sometable As g2   
WHERE g1.gid = 1 and g1.gid <> g2.gid  
ORDER BY distance(g1.the_geom,g2.the_geom) 
LIMIT 5

Now for moderate to large datasets, this exercise becomes painfully slow because distance is not an indexable function because it involves relations between two entities.

If you knew something about your data like there are no neighbors past 300 units away from your geometry of interest or that at that point you could care less what lies past 300 units, then you could significantly improve the speed of this query by writing it like this:

 SELECT g1.gid As gref_gid, g1.description As gref_description, g2.gid As gnn_gid, 
        g2.description As gnn_description  
    FROM sometable As g1, sometable As g2   
    WHERE g1.gid = 1 and g1.gid <> g2.gid AND expand(g1.the_geom, 300) && g2.the_geom  
    ORDER BY distance(g1.the_geom,g2.the_geom) 
    LIMIT 5

If you needed to run this for more than one record in sometable, then you can either put this in a loop and run it for each record of interest in sometable. Alternatively you could create an SQL function that returns a table for each record and use it like shown below.

CREATE OR REPLACE FUNCTION fnsometable_nn(gid integer, geom1 geometry, distguess double precision, n
umnn integer)  
    RETURNS SETOF sometable AS
    $BODY$
    SELECT g2.*  
        FROM sometable As g2   
            WHERE $1 <> g2.gid AND expand($2, $3) && g2.the_geom  
        ORDER BY distance(g2.the_geom, $2) LIMIT $4  
    $BODY$ LANGUAGE 'sql' VOLATILE;

To use the function to run for the first 100 records in sometable and return the 5 closest neighbors for each of those records, you would do this

SELECT g1.gid, g1.description, (fnsometable_nn(g1.gid, g1.the_geom, 300, 5)).gid As nn_gid, 
    (fnsometable_nn(g1.gid, g1.the_geom, 300, 5)).description as nn_description, 
    distance((fnsometable_nn(g1.gid, g1.the_geom, 300, 5)).the_geom,g1.the_geom) as dist
    FROM (SELECT * FROM sometable ORDER BY gid LIMIT 100) g1 

A lot of people are probably thinking "HUH" at this point. I will now wave my hand and leave it as an exercise for people to figure out why this works. I will say this much. I am using a feature in PostgreSQL which from my readings is an accident of design and only works with functions written in SQL and C language. This accident will probably be later corrected in future versions so don't expect this to work moving forward.

Returning tables in the select part of an sql statement is super non-standard and should I say bizarre and unelegant, but is the only way to do it currently in PostgreSQL where the you have a set returning function whose inputs depend on values from a relating table. This is because PostgreSQL currently does not support subselects and set returning functions in the FROM clause that dynamically change based on inputs from other table fields. In fact most databases do not support this.

Unfortunately PostgreSQL doesn't support a predicate similar to the CROSS APPLY predicate that SQL Server 2005 has. I think Oracle does by pretty much just leaving out the CROSS APPLY predicate. CROSS APPLY is more elegant looking I think in general and also more self-documenting in the sense that when you see the clause - you know what is being done.

SQL Server which I use very often - introduced this feature in the SQL Server 2005 offering and in a few cases its a life-saver.

I think the next next release of PostgreSQL will probably support a feature similar to this. Just to demonstrate how powerful this feature can be, if we had the cross apply feature, we could do away with our intermediate set returning function and simply write our query like this.

SELECT g1.gid, g1.description, nn.gid As nn_gid, nn.description As nn_description, distance(nn.the_g
eom,g1.the_geom) as dist
FROM (SELECT * FROM sometable ORDER BY gid LIMIT 100) g1 CROSS APPLY 
    (SELECT g2.*  
        FROM sometable As g2 
        WHERE g1.gid <> g2.gid AND expand(g1.the_geom, 300) && g2.the_geom)  
ORDER BY distance(g2.the_geom, g1.the_geom) LIMIT 5) As nn  

This would alleviate the need to create a new nn function for each table or kind of query we are formulating.

Its always interesting to see how different databases tackle the same issues. Oracle Locator and Oracle Spatial for example does have an SDO_NN and SDO_NN_DISTANCE functions which magically does all this for you, but I haven't used it so not sure of its speed or if it has any caveats. One I did notice from reading is that it doesn't deal with where clauses too well.



PostGIS Nearest Neighbor: A Generic Solution - Much Faster than Previous Solution

A generic solution to PostGIS nearest neighbor

After some heavy brainstorming, I have come up with a faster and more generic solution to calculating nearest neighbors than my previous solutions. For the gory details on how I arrived at this solution - please check out the following blog entries Boston GIS Blog entries on Nearest Neighbor Search.

In the sections that follow I will briefly describe the key parts of the solution, the strengths of the solution, the caveats of the solution, and some general use cases.

The parts of the solution

The solution is composed of a new PostgreSQL type, and 2 PostgreSQL functions. It involves a sequence of expand boxes that fan out from the bounding box of a geometry to as far out as you want. Below are the descriptions of each part.

  • pgis_nn - a new type that defines a neighbor solution - for a given near neighbor solution it gives you, the primary key of the neighbor and the distance it is away from the reference object
  • expandoverlap_metric - this is a helper function that returns the first fanned out box that an object is found in relative to the reference geometry. It returns an integer where the integer varies from 0 (in bounding box) to maxslices where maxslices is how many slices you want to dice up your larger expand box.
  • _pgis_fn_nn - this is the workhorse function that should never be directly called - it is a pgsql function that returns for each geometry k near neighbors
  • pgis_fn_nn(geom1 geometry, distguess double precision, numnn integer, maxslices integer, lookupset varchar(150), swhere varchar(5000), sgid2field varchar(100), sgeom2field varchar(100))) - this is our interface function written in language sql to fool the planner into thinking we are using an sql function. This is to get around certain anomalies of PostgreSQL that prevent pgsql set returning functions from being used if they reference fields in the from part of an SQL statement. Description of the inputs is as follows.
    1. geom1 - this is the reference geometry.
    2. distguess - this is the furthest distance you want to branch out before you want the query to give up. It is measured in units of the spatial ref system of your geometries.
    3. numnn - number of nearest neighbors to return.
    4. maxslices - this is the max number of slices you want to slice up your whole expand box. The more slices the thinner the slices, but the more iterations. I'm sure there is some statistically function one can run against ones dataset to figure out the optimal number of slices and distance guess, but I haven't quite figured out what that would be.
    5. lookupset - this is usually the name of the table you want to search for nearest neighbors. In theory you could throw in a fairly complicated subselect statement as long as you wrap it in parethesis. Something like '(SELECT g.*, b.name FROM sometable g INNER JOIN someothertable b ON g.id = b.id)' . Althought I haven't had a need to do that yet so haven't tested that feature.
    6. swhere - this is the where condition to apply to your dataset. If you want the whole dataset to be evaluated, then put 'true'. In theory you can have this correlated with your reference table by gluing parameters together, so can get pretty sophisticated.
    7. sgid2field - this is the name of the unique id field in the dataset you want to check for nearest neighbors.
    8. sgeom2field - this is the name of the geometry field in the dataset you want to check for nearest neighbors

The code can be downloaded here http://www.bostongis.com/downloads/pgis_nn.txt. I will also be posting this to the PostGIS Wiki once I clean it up a bit more.

Key features

My objective in creating this solution is to give PostGIS the same expressiveness that Oracle's SDO_NN and SDO_NN_DISTANCE functions provide Oracle. Unlike the Oracle solution, this solution is not implemented as an Operator, but instead as a set returning function.

Features are outlined below

  • Can handle more than one reference point at a time and return a single dataset. For example you can say give me the 5 nearest neighbors for each of my records in this query and get the result back as a single table.
  • You can limit your search to a particular distance range, slice it up as granularly as you want in order to achieve the best resulting speed. Note I implemented this function as linearly expanding boxes to utilize GIST indexes. In hindsight it might have been more efficient processor wise to have the boxes expand like a normal distribution or something like that.
  • You can apply where conditions to the sets you are searching instead of searching the whole table - e.g. only give me the 5 nearest hospitals that have more than 100 beds. You can also have the where conditions tied to fields of each record of interest. E.g. include all near neighbors, but don't include the reference point.
  • This solution works for all geometries - e.g. you can find point near neighbors to a polygon, polygon near neighbors to a polygon etc. Nearest distance to a line etc.

Caveats

  • This function will only work if you give it a primary key that is an integer. I figured an integer is what most people use for their primary keys and is the most efficient join wise.
  • This function currently doesn't handle ties. E.g. if you specify 5 neighbors and you have 6 where the last 2 are the same distance away, one will arbitrarily get left out. I think it will be fairly trivial to make it handle this case though.
  • If there are no neighbors for a particular reference that meet the stated criteria, that record would get left out. There are workarounds for this which I will try to provide in more advanced use cases.

Use Cases

All the data we will be using is in SRID: 26986. For a refresher course on loading data and setting up indexes, check Part 1: Getting Started With PostGIS: An almost Idiot's Guide

This returned 1000 records in 6 seconds on my dual windows 2003 server
Datasets used: buildings http://www.mass.gov/mgis/lidarbuildingfp2d.htm,
firestations: http://www.mass.gov/mgis/firestations.htm

Find nearest fire station for first 1000 buildings in Boston and sort in decreasing order of distance from the fire station such that the building furthest away from its nearest firestation is at the top.
SELECT g1.gid as gid_ref, f.name, g1.nn_dist/1609.344 as dist_miles, g1.nn_gid
FROM (SELECT b.gid, 
    (pgis_fn_nn(b.the_geom, 1000000, 1,1000,'fire_stations', 'true', 'gid', 'the_geom')).*  
    FROM (SELECT * FROM buildings limit 1000) b) As g1, fire_stations f
    WHERE f.gid = g1.nn_gid 
    ORDER BY g1.nn_dist DESC



Data set used: Fire stations from above and Massachusetts census tiger roads -http://www.mass.gov/mgis/cen2000_tiger.htm
For each street segment in zip 02109 list the 2 nearest fire stations and the distance in meters from each fire station. Order results by street name, street id and distance from nearest fire station
Total query runtime: 687 ms. 296 rows retrieved.

SELECT g1.gid as gid_ref, g1.street, g1.fraddl, g1.toaddl, f.name, g1.nn_dist as dist_meters, g1.nn_gid
FROM (SELECT b.*,
    (pgis_fn_nn(b.the_geom, 1000000, 2,1000,'fire_stations', 'true', 'gid', 'the_geom')).*  
    FROM (SELECT * FROM census2000tiger_arc where zipl = '02109') b) As g1, fire_stations f
    WHERE f.gid = g1.nn_gid 
    ORDER BY g1.street, g1.fraddl, g1.gid, g1.nn_dist 



For each street segment in zip 02109 list the 2 nearest buildings with building footprint > 5000 sq meters and specify how far away they are and the areas of each building.

--Total query runtime: 2656 ms. 296 rows retrieved.

SELECT g1.gid as gid_ref, g1.street, g1.fraddl, g1.toaddl, area(b.the_geom) as bldg_area, g1.nn_dist as dist_meters, g1.nn_gid
FROM (SELECT s.*,
    (pgis_fn_nn(s.the_geom, 1000000, 2,1000,'buildings', 'area(the_geom) > 5000', 'gid', 'the_geom')).*  
    FROM (SELECT * FROM census2000tiger_arc where zipl = '02109') s) As g1, buildings b
    WHERE b.gid = g1.nn_gid 
    ORDER BY g1.street, g1.fraddl, g1.gid, g1.nn_dist 


More Advanced Queries

Example of two simultaneous neighbor queries in one query

Layers used
http://www.mass.gov/mgis/schools.htm - schools
http://www.mass.gov/mgis/hospitals.htm - hospitals

Now for this, I think we can safely assume that no street we care about should be further than 30,000 meters from its nearest school and hospital so we will limit our search to 30,000 meters and dice up our box into 10 expand boxes - so boxes expand at a rate of 3000 meters out per expansion.

Give me the nearest elemenary school and the nearest hospital emergency center for street segments in select zipcodes and how far away they are from these street segments.
--Total query runtime: 5516 ms.(~5.5 secs). 291 rows retrieved.
SELECT emer.gid as gid_ref, emer.street, emer.fraddl, emer.toaddl, hospitals.name As emer_name, emer.nn_dist as emer_dist_meters, 
        schools.name as sch_name, sch.nn_dist as sch_dist_meters
FROM (SELECT b.gid, b.street,b.fraddl, b.toaddl,
    (pgis_fn_nn(b.the_geom, 30000, 1,10,'hospitals', 'er_status=''Y''', 'gid', 'the_geom')).*  
    FROM (SELECT * FROM census2000tiger_arc where zipl IN('02109', '02110')) As b) As emer INNER JOIN
    (SELECT b.gid,
    (pgis_fn_nn(b.the_geom, 30000, 1,10,'schools', 'true', 'gid', 'the_geom')).*  
    FROM (SELECT * FROM census2000tiger_arc where zipl IN('02109', '02110')) As b) As sch ON emer.gid = sch.gid
    INNER JOIN hospitals ON hospitals.gid = emer.nn_gid 
    INNER JOIN schools ON schools.gid = sch.nn_gid
    ORDER BY emer.street, emer.fraddl, emer.gid

Example limit distance exclusion queries

Provide a list of all schools in Massachusetts Boston extent whose closest hospital emergency room is further than 3000 meters away

Using our new neighbor function exclusion query
--Total query runtime: 2969 ms. 75 rows retrieved.
SELECT sch.gid as gid_ref, sch.name as sch_name
FROM schools sch LEFT JOIN (SELECT b.gid, 
    (pgis_fn_nn(b.the_geom, 3000, 1,2,'hospitals', 'er_status=''Y''', 'gid', 'the_geom')).*  
    FROM (SELECT * FROM schools WHERE schools.the_geom && GeomFromText('POLYGON((225446.34375 886435.5625,225446.34375 905270.8125,242278.640625 905270.8125,242278.640625 886435.5625,225446.34375 886435.5625))',26986)) b) As g1 ON sch.gid = g1.gid
    WHERE g1.gid IS NULL AND sch.the_geom && GeomFromText('POLYGON((225446.34375 886435.5625,225446.34375 905270.8125,242278.640625 905270.8125,242278.640625 886435.5625,225446.34375 886435.5625))',26986)
ORDER BY sch.name
The old way doing a simple expand distance exclusion query. SOMETIMES THE OLD WAYS ARE JUST BETTER
-- Total query runtime: 63 ms. 75 rows retrieved.
SELECT sch.gid as gid_ref, sch.name as sch_name
FROM schools sch LEFT JOIN 
(SELECT b.gid FROM schools b , hospitals h 
		WHERE b.the_geom && GeomFromText('POLYGON((225446.34375 886435.5625,225446.34375 905270.8125,242278.640625 905270.8125,242278.640625 886435.5625,225446.34375 886435.5625))',26986) AND h.er_status = 'Y'
				AND expand(b.the_geom, 3000) && h.the_geom AND distance(b.the_geom, h.the_geom) <= 3000) g1
ON sch.gid = g1.gid
WHERE g1.gid IS NULL AND sch.the_geom && GeomFromText('POLYGON((225446.34375 886435.5625,225446.34375 905270.8125,242278.640625 905270.8125,242278.640625 886435.5625,225446.34375 886435.5625))',26986)
ORDER BY sch.name

Caching data by storing closest neighbor in a table

For this exercise, we will create 2 new fields in our schools table called nearest_hospital and nearest_hospital_gid and update the fields with the closest hospital to that school. We will set our maximum distance search to 100,000 meters carved out into 10 10,000 meter expansions.

-- Query returned successfully: 2613 rows affected, 20625 ms (~1/3 minute) execution time.
ALTER TABLE schools ADD nearest_hospital character varying(255);
ALTER TABLE schools ADD nearest_hospital_gid integer;

UPDATE schools SET nearest_hospital = h.name , nearest_hospital_gid = h.gid
FROM (SELECT s.gid As sch_gid,
    (pgis_fn_nn(s.the_geom, 100000, 1,10,'hospitals', 'er_status=''Y''', 'gid', 'the_geom')).*  
    FROM schools As s) As nn_emer INNER JOIN hospitals h ON nn_emer.nn_gid = h.gid
WHERE schools.gid = nn_emer.sch_gid


Later on we realize we need an alternate hospital in case the first hospital is filled or because of traffic situation we can't get there, so we create a new set of fields to hold the second closest hospital. This particular exercise demonstrates where clause dependent on the additional attribute data of the reference record.


-- Query returned successfully: 2613 rows affected, 23032 ms (~1/3 minute) execution time.
ALTER TABLE schools ADD nearest_hospital_2 character varying(255);
ALTER TABLE schools ADD nearest_hospital_gid_2 integer;

UPDATE schools SET nearest_hospital_2 = h.name , nearest_hospital_gid_2 = h.gid FROM (SELECT s.gid As sch_gid, (pgis_fn_nn(s.the_geom, 100000, 1,10,'hospitals', 'er_status=''Y'' AND gid <> ' || CAST(s.nearest_hospital_gid As varchar(20)) , 'gid', 'the_geom')).* FROM schools As s) As nn_emer INNER JOIN hospitals h ON nn_emer.nn_gid = h.gid WHERE schools.gid = nn_emer.sch_gid AND nearest_hospital_gid IS NOT NULL


Part 1 - PostGIS and SharpMap in ASP.NET 2.0 using VB.NET: Compiling SharpMap with PostGIS

What is SharpMap

SharpMap is an opensource freely available mapping engine for the Microsoft.NET 2.0 Framework. It supports numerous datasources such as PostGIS, ESRI Shapefile, MSSQL Spatial, ECW and Oracle spatial. Future support is in the works for MapInfo Tab files and other datasources. Unfortunately as of yet it has not been thoroughly tested in Mono.NET and the SharpMap.UI (desktop portion) is noted to not work in Mono.Net - see notes for further details http://www.codeplex.com/Wiki/View.aspx?ProjectName=SharpMap&title=Can%20I%20use%20SharpMap%20on%20the%20Mono%20framework%20instead. What is especially interesting about SharpMap is that it is a relatively pure .NET implementation and can be used as both a Web mapping as well as a desktop toolkit.

If you are a .NET enthusiast as we are or just want to program in something like VB.NET, this is definitely something to take a look at.

Getting Started

In this exercise we will do what I call synchronized maps using PostGIS spatial database engine, SharpMap.Net mapping engine, and data from Boston Private Abandoned Property surveys. Synchronized Maps are maps laid out side by side that show the same location and are always in synch with each other. A common example of this is a keymap verses a full blown map where clicking in the key map zooms the main map and clicking on the main map zooms the key map. In our example we will be using maps that are all the same size, but represent different periods in time. I am going to use Visual Basic .NET for this example for a couple of reasons.

  • People rarely show VB.NET examples in .NET and focus on C#. This leaves VB.NET and general VBA developers like me feeling "What about me?"
  • C# doesn't have a With clause to my knowledge. How can people live without the convenience of a With clause :).
  • VB.NET reads more like a spoken language so I think easier to follow and debug.
  • VB.NET does a lot of transparent type-casting for you which C# doesn't. For example if you tried to glue together an int and a string, c# will yell (can't convert .. to ..) forcing you to write a mess of type-casting code where as VB.NET will observe that the integer can be converted to a string and do the casting automagically.

Pre-Requisites

  • Webserver or workstation with ASP.NET 2.0 installed
  • Microsoft.Net 2.0 Framework installed on workstation you will be using to compile
  • PostGreSQL 8.0 or above server with PostGIS 1.0 or above installed

Compiling the SharpMap Source

The default binaries of SharpMap available do not come precompiled with PostGIS support, so in this section, we'll go over how to compile your own

Before you can start, you'll need to gather the following items

  1. SharpMap.NET source - download the latest version from here, give it a .zip extension if it doesn't default to that when downloading, and unzip it unto your C drive. It should create a SharpMap folder.
  2. Copy the PostGIS.cs file from the extracted folder SharpMap/SharpMap.Extensions/Data/Providers to the folder SharpMap/SharpMap/Data/Providers.
  3. Download the latest PostgreSQL.NET driver (npgsql) binary for MS 2.0 from http://pgfoundry.org/projects/npgsql. For this I used Npgsql 1.0-bin-ms2.0.zip.

At this point, you can either use a development tool such as Visual Studio 2005, the freely available SharpDevelop ( http://www.sharpdevelop.com/OpenSource/SD/Default.aspx ) or Visual C# Express ( http://msdn.microsoft.com/vstudio/express/visualcsharp/ ) (a free download) and the included solution files to compile the source, or you can just compile with a command line using the .NET framework that you have already installed. I'm just going to describe how to do it with a commandline command since that doesn't assume any additional prerequisites and is simpler to explain.

  1. Create a bin folder in your extracted SharpMap/SharpMap folder and copy the Npgsql.* files and Mono.Security.dll from npgsql zip files into this new folder
  2. Creat a batch file in the bin folder with the following lines in it
    set croot=C:\SharpMap\SharpMap\
    "%SystemRoot%\microsoft.NET\Framework\v2.0.50727\csc.exe" /t:library /debug:full /out:%croot%bin\SharpMap.dll /recurse:%croot%*.cs /r:System.web.dll /r:System.data.dll /r:System.Xml.dll /r:system.dll /r:%croot%bin\Npgsql.dll
    pause

    Then run the batch script. Running the batch should create 2 files a .dll and a .pdb file. The pdb is used for debugging so it can highlight lines that break in the internal library.
  3. Now create an application folder on your webserver complete with a bin folder. Drop the SharpMap.dll, SharpMap.dbg, npgsql.dll, Mono.Security.dll files into the bin folder.

Loading the Test Data

I have packaged in the download file for this exercise, 2 sql scripts one to load neighborhoods, one for abandoned survey data. Below is the batch script to load all the data. For details on how to create your own load files, read Part 1: Getting Started With PostGIS: An almost Idiot's Guide or OGR2OGR Cheatsheet

Batch Script
c:\pgutils\psql -d gisdb -h localhost -U postgres -f neighborhoods.sql
c:\pgutils\psql -d gisdb -h localhost -U postgres -f abansurveys.sql
pause



Part 2 - PostGIS and SharpMap in ASP.NET 2.0 using VB.NET: Displaying the Maps

The Fun Part: Mapping the Data

There are 5 key things we need to do when presenting data on web maps. They are

  • Plotting the data
  • Maintaining view state - where was the user last - if this is there first visit - initialize the state.
  • Figuring out where a user clicked on an image map and converting these coordinates to spatial coordinates relative to view state.
  • Taking additional inputs such as request for specific layers or filtering of layers.
  • Doing something with these coordinates and inputs - e.g. do we Zoom in, Zoom Out, Pan, filter data


Everthing beyond those 5 basic items is presentational sugar.

Initializing Properties

For this exercise, we will use arrays to keep track of the 4 maps we have. This makes it easy to extend for more maps.

    Protected justmaps_string() As String = {"2002", "2003", "2004", "2005"}
    Protected justmaps_maps(3) As SharpMap.Map 'its a zeroth based array

The Presentation Layer

The key parts of our aspx is shown below. Note how we have named the image controls to be in line with our array variables. See the Demo in action

        <table>
            <tr>
                <td nowrap colspan="2">
                    <asp:RadioButtonList ID="rblMapTools" runat="server" RepeatDirection="Horizontal">
                        <asp:ListItem Value="0">Zoom in</asp:ListItem>
                        <asp:ListItem Value="1">Zoom out</asp:ListItem>
                        <asp:ListItem Value="2" Selected="True">Pan</asp:ListItem>
                    </asp:RadioButtonList>
                    <asp:Button ID="cmdReset" Text="Reset" runat="server"/>
                </td>
             </tr>
             <tr>
                <td nowrap>Unit Type: 
                    <asp:DropDownList ID="ddlLU" runat="server" AutoPostBack="true">
                        <asp:ListItem Value="ALL">ALL</asp:ListItem>
                        <asp:ListItem Value="A">A</asp:ListItem>
                        <asp:ListItem Value="C">C</asp:ListItem>
                        <asp:ListItem Value="CM">CM</asp:ListItem>
                        <asp:ListItem Value="I">I</asp:ListItem>
                        <asp:ListItem Value="R1">R1</asp:ListItem>
                        <asp:ListItem Value="R2">R2</asp:ListItem>
                        <asp:ListItem Value="R3">R3</asp:ListItem>
                        <asp:ListItem Value="RC">RC</asp:ListItem>
                    </asp:DropDownList>
                </td>
                <td>
                    Abandoned Type: 
                    <asp:DropDownList ID="ddlAbandtype" runat="server" AutoPostBack="true">
                        <asp:ListItem Value="ALL">ALL</asp:ListItem>
                        <asp:ListItem Value="A">Abandoned</asp:ListItem>
                        <asp:ListItem Value="B">Burned</asp:ListItem>
                        <asp:ListItem Value="D">Boarded</asp:ListItem>
                    </asp:DropDownList>
               </td>
            </tr>
            <tr>
                <td><b>Survey 2002</b><br />
                        <asp:ImageButton Width="300" Height="300" ID="imgMap2002" runat="server" style="border: 1px solid #000;" />
                </td>
                <td>
                    <b>Survey 2003</b><br />
                    <asp:ImageButton Width="300" Height="300" ID="imgMap2003" runat="server" style="border: 1px solid #000;" />
                </td>
            </tr>
            <tr>
                <td>
                    <b>Survey 2004</b><br />
                    <asp:ImageButton Width="300" Height="300" ID="imgMap2004" runat="server" style="border: 1px solid #000;" />
                </td>
                <td>
                    <b>Survey 2005</b><br />
                    <asp:ImageButton Width="300" Height="300" ID="imgMap2005" runat="server" style="border: 1px solid #000;" />
                </td>
            </tr>
        </table>

Initializing the Data

In ASP.NET the basic state is initialized in the Page Load event. Below is what our page load looks like

    Protected Sub Page_Load(ByVal sender As Object, ByVal e As System.EventArgs) Handles Me.Load
        Dim i As Integer
        Dim img As System.Web.UI.WebControls.Image

        For i = 0 To justmaps_string.GetUpperBound(0)
            img = Me.FindControl("imgMap" & justmaps_string(i))
            justmaps_maps(i) = Me.InitializeMap(New System.Drawing.Size(img.Width.Value, img.Height.Value), justmaps_string(i))
        Next

        If Page.IsPostBack Then
            'Page is post back. Restore center and zoom-values from viewstate
            For i = 0 To justmaps_string.GetUpperBound(0)
                justmaps_maps(i).Center = ViewState("mapCenter")
                justmaps_maps(i).Zoom = ViewState("mapZoom")
            Next
        Else
            'Save the current mapcenter and zoom in the viewstate by just picking the first map
            ViewState.Add("mapCenter", justmaps_maps(0).Center)
            ViewState.Add("mapZoom", justmaps_maps(0).Zoom)
            For i = 0 To justmaps_string.GetUpperBound(0)
                GenerateMap(justmaps_maps(i), justmaps_string(i))
            Next
        End If
    End Sub

Observe in the above, we create an instance of a sharpmap object to correspond with each .NET Image control. Below is what our initialize map code looks like.

    Public Function InitializeMap(ByVal size As System.Drawing.Size, ByVal maptype As String) As SharpMap.Map
        HttpContext.Current.Trace.Write("Initializing map...")
        'Initialize a new map of size 'imagesize'
        Dim map As SharpMap.Map = New SharpMap.Map(size)
        Dim ConnStr As String = ConfigurationManager.AppSettings("DSN")

        Dim dtItem As New SharpMap.Data.Providers.PostGIS(ConnStr, "abansurveys", "the_pointft")
        Dim dtNeighborhoods As New SharpMap.Data.Providers.PostGIS(ConnStr, "neighborhoods", "the_geom")
        Dim strWhere As String = "(yr = " & maptype & " AND abandtype <> 'N') "

        Dim lyrNeighborhoods As SharpMap.Layers.VectorLayer = New SharpMap.Layers.VectorLayer("Neighborhoods")
        Dim lyrNeighborhoodNames As SharpMap.Layers.LabelLayer = New SharpMap.Layers.LabelLayer("NeighborhoodNames")
        Dim lyrItem As SharpMap.Layers.VectorLayer = New SharpMap.Layers.VectorLayer("Item")

        If Not IsNumeric(maptype) Then 'map type should be an integer otherwise the request is a potential hack
            Return map
        End If
        With lyrNeighborhoods
            .DataSource = dtNeighborhoods
            .Style.Fill = Brushes.White
            .Style.Outline = System.Drawing.Pens.Black
            .Style.EnableOutline = True
        End With
        map.Layers.Add(lyrNeighborhoods)

        If Me.ddlLU.SelectedValue <> "ALL" Then
            strWhere &= "AND lu = '" & Me.ddlLU.SelectedValue & "' "
        End If

        If Me.ddlAbandtype.SelectedValue <> "ALL" Then
            strWhere &= "AND abandtype = '" & Me.ddlAbandtype.SelectedValue & "' "
        End If

        dtItem.DefinitionQuery = strWhere

        With lyrItem
            .DataSource = dtItem
            .Enabled = True
            With .Style
                .SymbolScale = 0.8F
            End With
        End With
        map.Layers.Add(lyrItem)

        With lyrNeighborhoodNames
            .DataSource = dtNeighborhoods
            .Enabled = True
            .LabelColumn = "name"
            .Style = New SharpMap.Styles.LabelStyle
            With .Style
                .ForeColor = Color.Black
                .Font = New Font(FontFamily.GenericMonospace, 11, FontStyle.Bold)
                .HorizontalAlignment = SharpMap.Styles.LabelStyle.HorizontalAlignmentEnum.Center
                .VerticalAlignment = SharpMap.Styles.LabelStyle.VerticalAlignmentEnum.Middle
                .Halo = New Pen(Color.Yellow, 2)
                .CollisionBuffer = New System.Drawing.SizeF(30, 30)
                .CollisionDetection = True
            End With
            .TextRenderingHint = Text.TextRenderingHint.AntiAlias
        End With
        map.Layers.Add(lyrNeighborhoodNames)

        map.BackColor = Color.White
        map.ZoomToExtents()
        HttpContext.Current.Trace.Write("Map initialized")
        Return map
    End Function

In the page load for each map we call a GenerateMap which renders the image of the map to the browser. The code looks like the below

    '<summary>
    'Grabs the map corresponding to a given image control, inserts it into the cache and sets the Image control Url to the cached image
    '</summary>
    Private Sub GenerateMap(ByVal aMap As SharpMap.Map, ByVal maptype As String)
        Dim img As System.Drawing.Image = aMap.GetMap()
        Dim imgID As String = SharpMap.Web.Caching.InsertIntoCache(1, img)
        CType(Me.FindControl("imgMap" & maptype), System.Web.UI.WebControls.Image).ImageUrl = "Getmap.aspx?ID=" + HttpUtility.UrlEncode(imgID)
    End Sub

Handling user requests

In this particular example, we've got a couple of actions that a user can perform.

  1. Zoom In
  2. Zoom Out
  3. Pan
  4. Filter by Unit Type
  5. Filter by Abandoned Type


To handle the first 3 actions, we define a single event handler that covers clicking on any of the maps. That code looks like

    Protected Sub imgMap_Click(ByVal sender As Object, ByVal e As System.Web.UI.ImageClickEventArgs) Handles imgMap2002.Click, imgMap2003.Click, imgMap2004.Click, imgMap2005.Click
        Dim i As Integer
        '-- Set center of the map to where the client clicked
        For i = 0 To justmaps_string.GetUpperBound(0)
            justmaps_maps(i).Center = justmaps_maps(i).ImageToWorld(New System.Drawing.Point(e.X, e.Y))
            '-- Set zoom value if any of the zoom tools were selected
            If rblMapTools.SelectedValue = "0" Then '//Zoom in
                justmaps_maps(i).Zoom = justmaps_maps(i).Zoom * 0.5
            ElseIf rblMapTools.SelectedValue = "1" Then '//Zoom out
                justmaps_maps(i).Zoom