Archive for the ‘Mapping’ Category

WMS GetFeatureInfo spanning the dateline in OpenLayers

Friday, July 27th, 2012

If you’re using OpenLayers to draw a WMS layer and calling GetFeatureInfo, you’re probably having trouble when looking across the antimeridian, where 180° longitude becomes -180°. This is because OpenLayers sends a bounding box to the WMS server (I’m only dealing with Geoserver 2.1.4 at the moment, not sure on other systems) that includes coordinates beyond the normal limit. If centered over New Zealand it might send a bounding box stretching from 160° to 190°. If your data is in normal (-180->180) space Geoserver won’t find any of the data beyond 180 and if you convert OpenLayers to sending rationalized coordinates (160° to -170°) Geoserver will complain that the minX is greater than the maxX. So what are we to do?

The solution I have come up with is to trim the box you are sending to the server so that it doesn’t cross the dateline. depending on where you were clicking in my example you’d either send 160° to 180° or -180° to -170°. The trick is to work out which side of the dateline you clicked on and which side your view is based on. If your view had been centered on the Chatham Islands OpenLayers would see those same boundaries as -200° to -170°.

Starting from the OpenGeo code to query a wms layer using OpenLayers I added come calculation to trim the bounding box and shift the X and Y values (which are the point the user clicked relative to the top left of the bounding box) if the left boundary has been moved.


map.events.register('click', map, function (e) {
bbox = map.getExtent();
width = map.size.w;
clickX = e.xy.x;

//Start dateline solution
clickRatio = e.xy.x/map.size.w;
fullMapWidth = Math.abs(map.getMaxExtent().left) + Math.abs(map.getMaxExtent().right);

if (bbox.right > map.getMaxExtent().right)
{
    datelineRatio = (map.getMaxExtent().right-bbox.left)/(bbox.right-bbox.left);

    if(datelineRatio < clickRatio)//Click was far side of dateline, normal bbox would not find data
    {
        bbox.left = map.getMaxExtent().left;
        bbox.right -= fullMapWidth;
        clickX -= width * datelineRatio;
        width = width * (1-datelineRatio);
    }

} else if (bbox.left < map.getMaxExtent().left){
    datelineRatio = (map.getMaxExtent().left-bbox.left)/(bbox.right-bbox.left);
    if(datelineRatio > clickRatio)//Click was far side of dateline, normal bbox would not find data
    {
        bbox.right = map.getMaxExtent().right;
        bbox.left += fullMapWidth;
    width = width * datelineRatio;
    }
}

//end dateline solution

var url = wfsPath
+ "?REQUEST=GetFeatureInfo"
+ "&EXCEPTIONS=application/vnd.ogc.se_xml"
+ "&BBOX=" + bbox.toBBOX()
+ "&X=" + Math.round(clickX)
+ "&Y=" + Math.round(e.xy.y)
+ "&INFO_FORMAT=text/html"
+ "&QUERY_LAYERS=" + wfsLayer
+ "&LAYERS=" + wfsLayer
+ "&FEATURE_COUNT=50"
+ "&SRS=EPSG:900913"
+ "&STYLES="
+ "&WIDTH=" + Math.round(width)
+ "&HEIGHT=" + map.size.h;

window.open(url,"getfeatureinfo","location=0,status=0,scrollbars=1,width=600,height=150");

});

Geotiffs: Stage 1 seamless USGS map tile sets

Tuesday, June 7th, 2011

For a little while now I’ve been manipulating the excellent USGS geotiffs to fed my PCT mapping fetish. Since I’ve worked out a good enough system I thought I might share it with the world.

First off you’ll need a bunch of stuff installed on your computer.

With those installed you’re ready to follow my process.

1. Remove the collar

The USGS maps that I’ve been using come with a ‘collar’, that bit around the outside of the actual map data that you’d have to cut off if you were printing them out and sticking them together.

USGS maps with collars

Before cutting the collar off make sure to keep a copy of the original tiff because most image editing programs will lose the geo data when you press Save. In Photoshop I had to convert the image to RGB mode so that I could make the outside part transparent. It’s important not to change the size of the image. The embedded coordinates are for the corner of the tiff even though the map doesn’t reach there. Save the image and get back to the command line for step 2.

USGS maps with collars

2. Restore the geospatial data

The new tiff you have just saved does not have the embedded coordinates that it started with. So long as you kept a copy of the original file you can reinstate the data using the following command and the gdalcopyproj program.

> gdalcopyproj.py original_file.tif new_file.tif

3. Create a Virtual Tile Set

If you’re combining multiple maps you’ll either need to merge them into one large file or create a virtual tile set like this

>gdalbuildvrt -srcnodata b4 -hidenodata merged.vrt *.tif

The “b4″ tells GDAL that band four of your images is the transparency layer, and the “hidenodata” tells it that when two maps overlap, hide the one that has no data (is transparent). Without that the transparent collars of the geotifs erase useful content from geotifs they join up with.

You can merge the files like this

> gdalwarp -co COMPRESS=LZW *.tif merged.tif

But it won’t do the compression until the end, so with USGS maps you’ll be generating files over 1GB even if you only have a few source maps. The virtual tile set method is definitely quicker and takes less disk space.

4. Cutting the tiles

The following command will take either your merged geotif file, or the virtual tile set as the first parameter and an output directory as the second (which will be created if it doesn’t exist).

>gdal2tiles.py merged.vrt tiles

The process can take quite a while, but it starts with the most detailed layer and each one after that is about 4 times faster than the previous one. When it’s complete you’ll have a full TMS tile stack ready for use on web map systems such as Google Maps and OpenLayers and my iOS app Tiled Maps. Next time I’ll post the scripts I use to optimise a tile set down from ~1.7Gb to 700Mb.