Friday 28 April 2017

postgis - Heat map / Density map from dynamic points table in MapServer (GeoServer)?


How can I create a heat map (some would call a density map) from a points layer in MapServer (or GeoServer)? (preferably able to choose Interpolation method, NN, IDW, etc and color map with a transparent color)


The points are stored in PostGIS table and the table is dynamic, meaning I need to create the heat map on the fly.



Answer



Mapserver can do it with the current development version which will be soon released as v7.0. How it works is best documented here: http://mapserver.org/development/rfc/ms-rfc-108.html


For testing the heatmap feature install MapServer 6.5-dev into your environment. Next download "pnts" shapefile (.shp, .shx, .dbf, and .prj) from github https://github.com/mapserver/msautotest/tree/master/gdal/data. Store the shapefile somewhere in your disk. Save the following mapfile as "heatmap.map" (slightly edited from https://github.com/mapserver/msautotest/blob/master/gdal/heat.map) into the same directory:



map
size 1000 500
extent -180 -90 180 90
name "test heat"
imagetype "png"
units dd
web
metadata
"ows_srs" "epsg:4326 epsg:3857 epsg:900913"
"ows_enable_request" "*"

end
end
projection
"+init=epsg:4326"
end
CONFIG "MS_ERRORFILE" "stderr"
layer
name "heatmap"
type raster
connectiontype kerneldensity

connection "points"
status on
processing "RANGE_COLORSPACE=%color%"
processing "KERNELDENSITY_RADIUS=%radius%"
processing "KERNELDENSITY_COMPUTE_BORDERS=%border%"
processing "KERNELDENSITY_NORMALIZATION=%norm%"
offsite 0 0 0
SCALETOKEN
NAME "%radius%"
VALUES

"0" "15"
"255000000" "20"
END
END
SCALETOKEN
NAME "%border%"
VALUES
"0" "ON"
"255000000" "OFF"
END

END
SCALETOKEN
NAME "%norm%"
VALUES
"0" "AUTO"
"255000000" "30"
END
END
SCALETOKEN
NAME "%color%"

VALUES
"0" "HSL"
"255000000" "RGB"
END
END
class
style
COLORRANGE "#0000ff00" "#0000ffff"
DATARANGE 0 32
end

style
COLORRANGE "#0000ffff" "#ff0000ff"
DATARANGE 32 255
end
end
end
symbol
name "circle"
type ellipse
points 1 1 end

end
layer
name "points"
status on
type POINT
data "pnts.shp"
projection
"+init=epsg:4326"
end
CLASS

MAXSCALE 255000000
STYLE
SIZE [VAL]
END
END
CLASS
MAXSCALE 265000000
STYLE
SIZE 0.1
END

END
CLASS
MAXSCALE 275000000
EXPRESSION ([VAL]>1)
STYLE
SIZE 1
END
END
CLASS
MAXSCALE 275000000

STYLE
SIZE 2
END
END
end
end

The following request, as edited to suit your installation, should show a heatmap on your browser:


http://localhost/cgi-bin/mapserv.exe?map=/ms4w/heatmap.map&mode=map&layers=heatmap


Heatmap is also available as WMS service at URL


http://localhost/cgi-bin/mapserv.exe?map=c:\ms4w\heatmap.map

Unfortunately I can't add a nice screen capture from my own computer here because I could not make it to work on Windows. Hopefully you have Linux and better luck. However, you can check some expected results from https://github.com/mapserver/msautotest/tree/master/gdal/expected, for example this https://github.com/mapserver/msautotest/blob/master/gdal/expected/heatmap-r20-noborder-fixednorm-rgb-expression.png


When you have managed to get the demo heatmap to work you will only need to edit the mapfile to read "points" layer from PostGIS instead of a shapefile and, if needed, to edit projection and extents to suit your data.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...