Storing Elevation data in GeoServer

We needed a way to store elevation data in GeoServer. Turns out this is very easy, just use Gdal to create greyscale floating point GeoTIFF files. You get all the Gdal raster manipulation, reprojection, resampling etc etc, while keeping all the Z accuracy in the floating point band. Then, you can use the Raster SLD to render as color ramps, or probably in the near future, hill shades.

In our case, we wanted to be able to handle very large areas of high precision Lidar data. I wrote an enhancement to our "Raster Chopper" that converts any Gdal supported input source to geoTiff grey-scale floating point quarter quad tiles in geographic 4326. Almost exactly like we are storing imagery. To display, we use the GeoTiff mosaicing GeoServer driver. Works nicely, and should scale to very large, very accurate Lidar data in the future.

Geoserver WMS Tiles in Silverlight Map Control

While Randy (see previous blog) showed how to use GeoWebCache to display tiles in the Silverlight map control, we needed to display WMS tiles that are not cached. The reason was that GWC doesn't play nice with GeoServer password secured layers. This would also be useful for layers that change frequently and cannot be cached.

I put together this TileSource class in VB.NET to do this:

'Normal WMS (non GWC)
'Translated from http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/
'
Public Class NormalWMSTileSource
Inherits Microsoft.VirtualEarth.MapControl.TileSource

Sub New()
MyBase.UriFormat = "http://ogi.state.ok.us/geoserver/wms?LAYERS=ogi:okcounties&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&EXCEPTIONS=application%2Fvnd.ogc.se_inimage&SRS=EPSG%3A900913&BBOX={0},{1},{2},{3}&WIDTH=256&HEIGHT=256&TRANSPARENT=TRUE"
End Sub

Sub New(ByVal WMSLayerName As String)
MyBase.UriFormat = "http://ogi.state.ok.us/geoserver/wms?LAYERS=" & WMSLayerName & "&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&EXCEPTIONS=application%2Fvnd.ogc.se_inimage&SRS=EPSG%3A900913&BBOX={0},{1},{2},{3}&WIDTH=256&HEIGHT=256&TRANSPARENT=TRUE"
End Sub


Public Overrides Function GetUri(ByVal x As Integer, ByVal y As Integer, ByVal zoomLevel As Integer) As System.Uri
'Example:
'http://ogi.state.ok.us/geoserver/wms?LAYERS=ogi%3Aokcounties&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&EXCEPTIONS=application%2Fvnd.ogc.se_inimage&SRS=EPSG%3A900913&BBOX=-11271098.44125,4070118.8815625,-10958012.373437,4383204.949375&WIDTH=256&HEIGHT=256


'Convert x,y,zoom to 900913 coords and return a wms URI
Dim R As Rect = TileBounds(x, y, zoomLevel)

Return New Uri(String.Format(Me.UriFormat, R.Left, R.Top, R.Right, R.Bottom))
End Function


'def(__init__(self, tileSize = 256))
'"Initialize the TMS Global Mercator pyramid"
'self.tileSize = tileSize
'self.initialResolution = 2 * Math.PI * 6378137 / self.tileSize
'# 156543.03392804062 for tileSize 256 pixels
'self.originShift = 2 * Math.PI * 6378137 / 2.0
'# 20037508.342789244


Dim tileSize = 256
'"Initialize the TMS Global Mercator pyramid"
Dim initialResolution = 2 * Math.PI * 6378137 / tileSize
'# 156543.03392804062 for tileSize 256 pixels
Dim originShift = 2 * Math.PI * 6378137 / 2.0
'# 20037508.342789244


'def(TileBounds(self, tx, ty, zoom))
'"Returns bounds of the given tile in EPSG:900913 coordinates"

'minx, miny = self.PixelsToMeters(tx * self.tileSize, ty * self.tileSize, zoom)
'maxx, maxy = self.PixelsToMeters((tx + 1) * self.tileSize, (ty + 1) * self.tileSize, zoom)
'return ( minx, miny, maxx, maxy )

Private Function TileBounds(ByVal tx, ByVal ty, ByVal zoom) As Rect

'"Returns bounds of the given tile in EPSG:900913 coordinates"

Dim MinP As Point = PixelsToMeters(tx * tileSize, ty * tileSize, zoom)
Dim MaxP As Point = PixelsToMeters((tx + 1) * tileSize, (ty + 1) * tileSize, zoom)

Dim R As New Rect(MinP, MaxP)

Return R

End Function

'def(PixelsToMeters(self, px, py, zoom))
'"Converts pixel coordinates in given zoom level of pyramid to EPSG:900913"

'res = self.Resolution(zoom)
'mx = px * res - self.originShift
'my = py * res - self.originShift
'return mx, my

Private Function PixelsToMeters(ByVal px, ByVal py, ByVal zoom) As Point

Dim P As New Point

Dim res = Resolution(zoom)
P.X = px * res - originShift

'Strange, this come out negative? Should be positive.
P.Y = -(py * res - originShift)

Return P
End Function

'def(Resolution(self, zoom))
'"Resolution (meters/pixel) for given zoom level (measured at Equator)"

'# return (2 * math.pi * 6378137) / (self.tileSize * 2**zoom)
'return self.initialResolution / (2**zoom)

Private Function Resolution(ByVal zoom)
'"Resolution (meters/pixel) for given zoom level (measured at Equator)"

'# return (2 * math.pi * 6378137) / (self.tileSize * 2**zoom)
Return initialResolution / (Math.Pow(2, zoom))

End Function



End Class

Black areas in Silverlight Map Control using GeoWebCache

Strange behavior in Silverlight map control (SMC).

This blog from Randy George shows how to use GeoServer GeoWebcached layers in the Silverlight map control:

http://www.cadmaps.com/gisblog/?p=54

Very nice, works well. I used this syntax for my Tile Source

MyBase.UriFormat = http://ogi.state.ok.us:8080/geoserver/gwc/service/gmaps?layers=ogi:okcounties&zoom={2}&x={0}&y={1}

which pulls tiles from the GeoWebCache.

Problem is, some of our layers have a smaller coverage area than others. GeoWebCache is smart enough to not try to generate tiles if the layer doesn't cover the area. Problem is, that the SMC shows a tile that doesn't return data as BLACK. Very ugly - see picture.

So, how to fix? Likely this is a bug in the SMC, so it needs to be fixed there. But in the meantime, you can change the extents of your GeoServer layer to be larger. Then, reload the GeoWebCache by going to the http://yourserver:8080/geoserver/gwc/demo and clicking on "Reload Configuration" - or possibly restarting GeoServer, I didn't try that.

Then GWC will send back transparent PNGs instead of nothing at all, and the SMC is happy.

Adding Personalization to ASP.NET with PostGIS (PostgreSql)

Next step - adding personalization to the web site using Postgresql.

Very easy, as well as documented.

Just go to:

http://dev.nauck-it.de/aspsqlprovider/

and follow the directions.

One thing not mentioned is the PostgreSQL encoding of the database when you create it. I used the default Win1252, which seemed to work fine.

I created the database, ran the SQL script to create the tables, and added/modified the web.config file.

Created a asp.net web form with a login control, ran it, but it didn't work. Got:

at NauckIT.PostgreSQLProvider.PgMembershipProvider.ValidateUser(String username, String password) in C:\Users\Daniel Nauck\Documents\Visual Studio 2008\Projects\PostgreSQLProvider\PgMembershipProvider.cs:line 1229 at System.Web.UI.WebControls.Login.AuthenticateUsingMembershipProvider(AuthenticateEventArgs e)

Well, that was a bit disconcerting. So, looked at the Immediate window in VS2008, for more info, and saw:

Npgsql.NpgsqlException:
syntax error at or near "("
Severity: ERROR
Code: 42601
at Npgsql.NpgsqlState.d__a.MoveNext()
at Npgsql.NpgsqlState.IterateThroughAllResponses(IEnumerable`1 ienum)
at Npgsql.NpgsqlState.ProcessBackendResponses(NpgsqlConnector context)
at Npgsql.NpgsqlReadyState.Flush(NpgsqlConnector context)
at Npgsql.NpgsqlConnector.Flush()
at Npgsql.NpgsqlCommand.Prepare()
at NauckIT.PostgreSQLProvider.PgMembershipProvider.ValidateUser(String username, String password) in C:\Users\Daniel Nauck\Documents\Visual Studio 2008\Projects\PostgreSQLProvider\PgMembershipProvider.cs:line 1210
A first chance exception of type 'System.Configuration.Provider.ProviderException' occurred in NauckIT.PostgreSQLProvider.DLL

Not too helpful...

Then tried various things, made sure the uid and password was ok, changed the name of the application to the name of my project, etc. These had no effect.

Then, did a get of the latest Nauck from SVN and added to my solution. Then it started working. Perhaps the release 1.240 is a little bogus? Anyway, now it works. Awesome.

Converting a Lat Long to a USGS Quarter Quad Number in VB

USGS uses an odd system of numbering the 24k quads, and quarter quads.
it uses lat, lon(lower RIGHT corner), then letters A - H from right to left, and 1 -8 from bottom to top. Quarter quads are then numbered 1 is upper right, 2 is upper left, 3 is lower left, 4 is lower right of the whole quad. eg 33106-E32 where the 2 at the end is the qq number.

An irritating little problem, but here is a solution in VB.net


Given a lat-lon, what is the USGS quad number (or USGS quarter quad in this case)

Private Function LatLonToQQName(ByVal Latitude As Double, ByVal Longitude As Double) As String

'What is the Longitude whole number part (degree) of the right corner
Dim LonDegree As Integer = Math.Ceiling(Longitude)

'Get the fractional part of lon
Dim FLon = LonDegree - Longitude

'Calculate the Number, from 1 to 8, starting from right corner
Dim Number As String = Math.Ceiling(FLon / 0.125).ToString

'Next, the number, from 1 to 8, starting from bottom
Dim LatDegree As Integer = Math.Floor(Latitude)

'Get the fractional part of lon
Dim FLat = Latitude - LatDegree

Dim Letter = Chr(Asc("A") + Math.Floor(FLat / 0.125)).ToString

'Finally, the QQ number (from 1 to 4) 1 is upper right, 2 is upper left, 3 is lower left, 4 is lower right
Dim QQNumberLon As Integer = Math.Floor((FLon / 0.125 - Math.Floor(FLon / 0.125)) * 2)
Dim QQNumberLat As Integer = Math.Floor((FLat / 0.125 - Math.Floor(FLat / 0.125)) * 2)

Dim QQNumber As String = ""
If QQNumberLat = 0 And QQNumberLon = 0 Then
QQNumber = "4"
ElseIf QQNumberLat = 0 And QQNumberLon = 1 Then
QQNumber = "3"
ElseIf QQNumberLat = 1 And QQNumberLon = 1 Then
QQNumber = "2"
ElseIf QQNumberLat = 1 And QQNumberLon = 0 Then
QQNumber = "1"
End If

Return LatDegree.ToString & (Math.Abs(LonDegree)).ToString("D3") & "-" & Letter & Number & QQNumber
End Function