Azavea Labs

Where software engineering meets GIS.

Restricting Zoom with Multiple OL Basemaps

DistrictBuilder logoAs David recently posted, our team has been hard at work implementing DistrictBuilder, where we’ve been investing a great deal of effort on both performance and usability. One feature we added in the spring was the addition of basemaps to the user interface. Before this addition, users labored over drawing the perfect district configurations without a whole lot of context of the surrounding environment (e.g. roads, water boundaries, etc.). When the time came to add a basemap to the application, it didn’t feel right restricting it to a single type of map, or even a single provider. We wanted to allow for users to have the choice to select the best map for the task at hand. Could an application promoting democracy really have it any other way?

We set out to support several base map options as well as any combination of options, including:

  • Bing Maps (satellite, roads and hybrid)
  • GoogleMaps (satellite, roads and hybrid)
  • ArcGIS Online (any of several maps)
  • OpenStreetMap

Since DistrictBuilder needs to be flexible enough to meet the needs of users and administrators in a variety of situations, we decided on a two step approach to basemap configuration. First, the administrator specifies, in the configuration file, which of the combinations of map providers and map types are allowed to be selected. Then DistrictBuilder presents all of the configured options to the user, where they can be toggled among at any time while a plan is being viewed or edited.

Here’s an example of an instance configured with an OpenStreetMap road layer, a Bing hybrid layer, and a Google satellite layer:

Road, Hybrid, and Satellite

Here’s another example with only road layers — one for each of the three configured providers:

Roads for three vendors

DistrictBuilder currently allows the configuration of basemaps using permutations of each of the three vendors and three map types described above. Adding more options is a relatively easy task, however. With the launch of Fix Philly Districts, we wanted the basemap colors to be slightly more muted than the above options, and ended up adding support for the ArcGIS Online World Topographic Map. We also experimented with the Google Maps V3 custom styling API, which looked great, but introduced performance problems when panning and zooming (animations).

There were, of course, some hoops that needed to be jumped through in order to get all of these basemaps behaving correctly on the same map, which will be discussed below. I’ve extracted the logic required to do so into a small demo that can be viewed/downloaded here. The demo has also been embedded into this post, and can be interacted with without going anywhere:

Zoom Levels

Many of the challenges that needed to be overcome to get this working correctly were brought about because we needed to restrict the zoom levels to the area at hand. We wanted to eliminate superfluous zoom levels to ensure the user was always operating within the appropriate boundaries (note: it is not done in this demo, but in DistrictBuilder we also restrict the extent with the ‘restrictedExtent’ map parameter, so users can’t even pan outside of the area).

One difficulty with setting zoom levels on the different layers is that the layers don’t use zoom parameters consistently. In Bing (the VirtualEarth layer), minZoomLevel and maxZoomLevel are needed. In Google, minZoomLevel is needed, but it requires numZoomLevels instead of maxZoomLevel. And in OpenStreetMap (the OSM layer), well…no combination of those seem to work — we needed to slightly modify the XYZ layer (OSM’s base class) to allow maxResolution to be changed based on the minZoomLevel. To see how this is done, view the demo source. With that change in place, the list of required layer parameters is as follows:

  • Bing – minZoomLevel, maxZoomLevel, projection, sphericalMercator, maxExtent
  • Google – numZoomLevel, minZoomLevel, projection, sphericalMercator, maxExtent
  • OpenStreetMap - numZoomLevel, minZoomLevel, projection

Coordinate Systems

We also faced some problems related to coordinate systems. DistrictBuilder uses GeoServer and GeoWebCache to serve up WMS layers. The coordinate system of our data is one version of the the ever-changing “Popular Visualization CRS / Mercator” projection. We needed to match up the OpenLayers projection to the one used on our data, or else we were seeing slight offsets on our overlays. Unfortunately, the ‘projection’ layer parameter isn’t always used within the layers correctly. For example, any layer using the SphericalMercator class gets its projection automatically hardcoded to 900913. We needed to make a slight modification to the SphericalMercator class to allow the ‘projection’ parameter to carry through. This can be seen by viewing the demo source.

Bonus: Math Time!

One interesting part about implementing zoom restriction was that we needed it to work in any instance of DistrictBuilder — from large states to small towns, which may have vastly different extents. Instead of having an administrator figure out the proper minimum zoom level, we calculate it automatically based on the extent, which requires a little bit of basic algebra.

For Philadelphia, the extent of our area is:

[-8397913.926216, 4842467.609439, -8329120.600772, 4895973.529229]

In DistrictBuilder, we calculate this dynamically on the server side (using Django) by filtering all of the geounits in the database and calling the ‘extent’ function on the query set. For the demo, this is hardcoded. Here’s how to transform this extent into a Spherical Mercator minZoomLevel:

  • Find the width of the area in meters.
var studyWidthMeters = extent[2] - extent[0];
  • Find the width of the map in pixels. In the demo, this is hardcoded, because we are setting the div size of the map. In DistrictBuilder, the map takes up the whole screen, and this value is calculated on the fly based on the size of the div in which the map occupies.
var mapWidthPixels = 450;
  • Find the map resolution, or meters per pixel.
var resolution = studyWidthMeters / mapWidthPixels;
  • Find the maximum map resolution. In Spherical Mercator, the maximum resolution is one 256×256 tile taking up the entire circumference Earth. So dividing the circumference of the earth (~40,000km) by 256 gives us the maximum meters per pixel, which is a constant.
var maxResolution = 156543.033928;
  • Spherical Mercator zoom levels work like a pyramid. Each zoom breaks the current tile up into a 2×2 group of 256×256 tiles, essentially halving the resolution each time. Therefore, finding the resolution at a given zoom level looks like this:
maxresolution / 2^zoom = resolution
  • We know the resolution and max resolution already, and need to find the zoom:
zoom = log(maxresolution/resolution)/log(2)
  • Or in javascript:
var minZoom = Math.log(maxResolution / resolution) / Math.LN2;

Building Districts in Web-Time

DistrictBuilder logoMost recently, the Politics, Redistricting and Elections team has been working closely with the Public Mapping Project to build DistrictBuilder, an open source, web-based application that enables regular citizens to use powerful tools to draw their own legislative districts. If you’ve seen how badly the professionals can mangle districts (Exhibit A, Exhibit B, etc), it’s easy to imagine that any given citizen, given the right tools, could do it better.

We spent quite a bit of time making the application easy to use and responsive in modern desktop web browsers.  The “easy to use” part was tackled by our excellent UI/UX design team. The “responsive” part was the domain of  our engineers.  That’s where the fun began for me.

DistrictBuilder is designed to use any polygon shapefile, transform it into an internal data model, then make that accessible via map tiles and geometric features.  When serving map tiles, we use GeoServer and GeoWebCache to generate the tiles and cache them, respectfully. This performance is great — pre-generated map tiles are the best we can aim for with respect to the base map tiles. Serving geometric features at full resolution, however, introduces a slew of problems. A few that stood out right away:

  • Web Browser Limitations — 9 out of 10 experts agree: too many map features has a significant performance impact on web browsers, with the greatest impact on the Microsoft Internet Explorer browser.
  • Excessive Coordinates — delivering lots of polygon coordinate pairs that the user may never see consumes valuable bandwidth and rendering time.
  • Server Processing Time — recalculating state-wide geometric features consumes valuable CPU time.

Web Browser Limitations

First, we tackled the browser performance issues. A sluggish browser is the kiss of death in the web world, and we had to make the application experience as fast as possible before looking at the server processing time.

We originally gave users the power to create highly detailed districts at the statewide level, but realized that no modern web browser could handle the volume of polygon features that would need to be served to represent an entire state.  In order to mitigate this limitation, we limited the size and number of features sent to the browser. With some scale-dependent logic, a user zoomed in to a detail of a district can finely tune the boundary by moving smaller geographic features (e.g. census blocks), and a user zoomed out to the state-wide level can manipulate the districts by moving large geographic features only (e.g. counties). In addition, when editing the finest details, we limit the number of features a user can move in a single edit.

Excessive Coordinates

The next thing to go was the set of full resolution geometries. In DistrictBuilder, users never actually see the full geometries, but an adaptively simplified (sometimes called generalized) geometry; depending on the scale of the map view, the server will deliver geometries with appropriate coordinate resolutions. Simply put: as you zoom in on the map, you get more detail in the geometries.

By simplifying counties, the geometries are reduced from 166,958 points to 4,821. When a user is zoomed out, there is no noticeable difference between these geometries!  However, as the user is interacting with higher resolution maps, DistrictBuilder loads in higher-resolution geometries on demand. The following images demonstrate the difference in the geometry detail:

Low Resolution Transition

The zoomed in County layer, with a low resolution district overlay (orange line). There are currently 1,414 coordinates in this view of the district overlay.

High Resolution Transition

The zoomed in VTD layer, with a high resolution district overlay (orange line). There are currently 3,253 coordinates in this view of the district overlay.

You can notice the differences in the district detail if you look closely at the orange district boundary. This transition happens seamlessly in the application, loading in the higher resolution geometries as web users zoom in to areas of interest.

We also eliminated coordinates that you never see.  It made no sense to serve  coordinates that were located in the opposite side of the state where a user was editing, just like you wouldn’t expect to get an encyclopedia in the mail when releasing an RFI. With the OpenLayers library, Strategies came in handy here, particularly BBOX.

Server Processing Time

After we had optimized the performance of the user interface, we shifted our focus to the server-side processing.  One of the features that makes DistrictBuilder such a powerful tool is the accuracy of the underlying data and constant feedback of important district statistics. In order to calculate all these statistics on the fly, it is necessary to leverage some tricks already mentioned with respect to map tiles: caching and generalizing.

Computation of the district statistics must happen every time a district boundary is changed. A naive solution to this problem would be to aggregate the values within the boundary every time a change is made.  This approach results in horrible performance. Instead, we just determine what has changed — which areas were added, which areas were removed — and recompute the delta, or change, on the previous district value.

Another trick to optimizing performance is in the way we determine the changing boundaries.  I’ll describe the problem using the census geographies of counties, tracts, and blocks. The structure and detail of the underlying data yielded computationally expensive queries against the block geometries.  We came up with a method of searching for the geographies in a hierarchical fashion — searching the counties first, then continuing to the next smallest-scale geography only if there was any remaining geometry left in the query.  We did the same for the tracts, and took a shortcut at the block level to exclude the block geometries.  This increased server side performance considerably.

King William County

King William county is comprised of 22 Voter Tabulation Districts and 1,527 Census Blocks.

Consider the following scenario: a user wants to move King William County (highlighted in yellow) from District 1, which is over populated, to District 3, which is under populated. Changing the boundaries with all the blocks in King William County would require testing at least 4,000 blocks for spatial intersections, then aggregating 1,527 data values, and recomputing the spatial aggregate (union) of those 1,527 geometries. With our hierarchical approach, we can change the boundary of the district with the county boundary, and change the population totals by the county’s population. A few orders of magnitude fewer operations to perform, and much faster from the user’s perspective.

Lessons Learned

Throughout the DistrictBuilder development process, the same core performance challenge has arisen: the volume of data must be reduced. This applies to all aspects of the application:

  • Map Tiles: pre-render tiles to keep the number of rendered tiles to a minimum at runtime.
  • Map Features: deliver to the browser only as much information as you can see (perhaps even less).
  • Database Queries: do anything possible to ensure that geometric operations are performed on simplified geometries.
  • Aggregating Statistics: cache whatever you can, and only compute the difference from the last cache state.

The above steps reduced the sheer number of operations and volume of processing that both the server and browser need to complete when creating new districts. These are lessons that translate well to any “big data” problem, and are crucial in bringing sophisticated GIS operations to the web.

Using the CQL_FILTER parameter with OpenLayers WMS layers

I’ve used Openlayer’s Marker layer in several projects and have always just accepted that I can’t display more than around 500 markers at a time for a given query. Recently, I found another way. We’re using GeoServer as a WMS tile server for the tree and municipal boundary layers in PhillyTreeMap.org. GeoServer’s WMS implementation allows an additional parameter in the url called CQL_FILTER. This parameter allows you to use a little language called Common Query Language, or CQL, to apply data filters to the tiles that GeoServer generates. CQL is a plain text, human readable query language created by the OGC, but I like to think of it as an extremely limited third-cousin-by-adoption of SQL. I haven’t found too much in the way of documentation on this obscure little gem, so here’s a rundown of how we used it to display search results in PhillyTreeMap.org.

If you look through the CQL and ECQL page in GeoServer’s documentation, there are several examples but they don’t cover everything you can do with CQL. Basically, a CQL filter is a list of phrases similar to where clauses in SQL, separated by the combiner words AND or OR. You can use the following operators in a CQL phrase:

  • Comparison operations: =, <, >, and combinations
  • Math expressions: +, -, *, /
  • NOT
  • IS, EXISTS
  • BETWEEN, BEFORE, AFTER
  • IN
  • LIKE
  • Geometric operators: CONTAINS, BBOX, DWITHIN, CROSS(ES), INTERSECT(S)

Some of these operations have examples in the GeoServer documentation, and others can be inferred from the GeoTools documentation (the stuff in their CQL.toFilter() calls). CQL can also call any of GeoServer’s filter functions.You can add parenthesis as needed to affect the order the filters are evaluated in, just like in SQL.

CQL has a lot of power for such a short spec, but it has a one very large deficiency that requires some database designing to avoid: the utter lack of join support. This makes sense when you consider that GeoServer doesn’t know about joins either. Ultimately, you’re using CQL against the GeoServer layer, not the underlying database structure. Building views or adding reference columns to the table GeoServer is accessing can help get around this.

In PhillyTreeMap.org, we use 4 types of CQL filters: id lookups using =, null checks using IS, date and integer ranges using BETWEEN and text searches using LIKE. Here are examples of those uses along with some array joining to get a valid CQL filter string at the end:

filter_list = []
filter_list.append("species_id = 212")
filter_list.append("height IS NULL")
filter_list.append("dbh BETWEEN 10 and 20")
filter_list.append("neighborhood_id_list LIKE '%42%' ")
cql = ' AND '.join(filter_list)
# should look like this: "species_id = 212 AND height IS NULL AND dbh BETWEEN 10 and 20 AND neighborhood_id_list LIKE '%42%' "

The above CQL filter would locate trees in a specific species, have no height value, only have dbh values between 10 and 20 inches and are in a particular neighborhood. The neighborhood_id_list filter would have been a join if this were written in standard SQL since neighborhoods and trees have a many-to-many relationship in our database. Since we can’t do joins, any time a tree is added or the location is updated, all of it’s related geographies’ ids are added to a reference column on the tree, and used specifically for this type of query.

CQL is passed to GeoServer in the same way as any other WMS variable. We’re using openlayers, so most of the WMS configuration variables are already set when we create the layer. The WMS layer has a little method called mergeNewParams that lets us change those parameters after the layer has been initialized. It also automatically redraws the layer, so the changes take place immediately. To add CQL to the WMS call, just add the CQL_FILTER variable to the layer’s parameters and the layer should update.

wms_layer.mergeNewParams({'CQL_FILTER': "species_id = 212 AND height IS NULL AND dbh BETWEEN 10 and 20 AND neighborhood_id_list LIKE '%42%' "});

You can remove any filters by deleting the parameter from the layer as if it was a normal javascript object. You’ll need to redraw the layer yourself before the change will be visible.

delete wms_layer.params.CQL_FILTER;
wms_layer.redraw();

jQuery, OpenLayers.Layer.Vector and IE8

One of the very first things you learn about OpenLayers is the importance if initializing maps in either the body tag’s onload event or in script tags at the bottom of the page. The basic reason for this is that many of the OpenLayers components need a page’s dom to be complete before it starts adding it’s own dom structure. The onload event is thrown after the browser is finished loading everything else, and the bottom of the page naturally gets rendered last, so these two places tend to work the best. Makes sense right?

If you’ve been using additional frameworks, like  jQuery ( or EXT or Django or Drupal or whatever) you’ve learned to rely on that framework’s onready event, or its rough equivalent.  This event tends to be thrown before the onload event so that the frameworks don’t need to wait for every last image and tag to be loaded. In theory, this could be a very problematic place to put OpenLayers map initialization code. The dom elements that OpenLayers is looking for may or may not be there, and there’s no way to tell. In practice, most browsers, and most OpenLayers components, work fine when initialized in a framework’s onready event. One layer in particular (OpenLayers.Layer.Vector) will not play nice with one particular browser (IE8) if initialized in a framework’s onready event. The vector layer depends on the document object’s namespace attribute, and this attribute simply isn’t always available before IE’s onload event (which is called after the onready event). As you can imagine, this was one headache I wish I could have avoided.

Thankfully, there is a fairly simple fix. Instead of the usual convention of :

$.ready(function() { init_my_map(); });

Add a stand-alone script tag to the extreme end of the body content on the pages where you need a vector layer:

<script>
    init_my_map();
</script>

OpenLayers meet ArcGIS.com

“Hey guys, we can support ArcGIS online basemaps, right?”
“Oh yeah, definitely, we use OpenLayers, there must be a layer for that”

Azavea is no stranger to contributing to open source projects, so when we learned that OpenLayers was lacking full support for the newer ArcGIS basemaps, I knew I had a job to do. The community had been looking for this functionality for over 2 years, and still didn’t have a solution. You might say that developing tools for crime fighters would be exciting and rewarding, but it is also nice when something you wrote has the potential to be used everyday by a community you value.

We quickly discovered that a few partial solutions existed, but something was ever so slightly off, those layers weren’t performing in all the configurations we needed. We had a head scratching situation, we had tiles, we had a standard format, but it seemed like they wouldn’t let us get the calibration just right. Finally the eureka moment hit when we discovered that the server was not being entirely literal with it’s capabilities, and that they were actually subtly changing as the zoom level changed! Once we knew that, it was a simple matter of overriding some base code, recalibrating and extrapolating the extent, and we were set!

Fast forward to today, and we’ve peer-reviewed the layer internally, added documentation, multiple examples, and unit tests, and we submitted a patch to the community about a month ago. Although it isn’t officially in the trunk yet, we already have word from a new friend in Spain that the layer has helped him with a project!

I’m excited that this layer could help more developers and organizations take advantage of the great cartography offered by Esri while using the powerful OpenLayers platform. So please keep your fingers crossed while we go through the approval process once again, or checkout the patch and try it for yourself, and let us know what you think. It can support out of the box ArcGIS.com basemaps, local ArcGIS servers, and raw file caches produced by ArcGIS.

Thank you to those who helped (Jeff, Aaron, Zwarg, and Robert), and to those who helped online and in the community, Thanks! Also thanks to Esri, who just made their ArcGIS map tiles available for free to everyone!