code

Saturday, 16 September 2017

MapServer OpenLayers Example.

The following is a simple example of using OpenLayers with a MapServer WMS.

It is based on the OpenLayers examples.

<!doctype html>
<html lang="en">
  <head>
    <link rel="stylesheet" href="ol.css" type="text/css">
    <style>
      .map {
        height: 600px;
        width: 100%;
      }
    </style>
    <script src="ol.js" type="text/javascript"></script>
    <title>OpenLayers (4.3.2) MapServer WMS - Wales Imagery Example</title>
  </head>
  <body>
    <h2>OpenLayers (4.3.2) MapServer WMS - Wales Imagery Example</h2>
    <div id="map" class="map"></div>
<div id="mouse-position"></div>
    <form>
      <label>Projection </label>
      <select id="projection">
        <option value="EPSG:4326">EPSG:4326</option>
      </select>
      <label>Precision </label>
      <input id="precision" type="number" min="0" max="12" value="4"/>
    </form>
    <script type="text/javascript">
      var mousePositionControl = new ol.control.MousePosition({
        coordinateFormat: ol.coordinate.createStringXY(4),
        projection: 'EPSG:4326',
        // comment the following two lines to have the mouse position
        // be placed within the map.
        className: 'custom-mouse-position',
        target: document.getElementById('mouse-position'),
        undefinedHTML: '&nbsp;'
      });
  
      var map = new ol.Map({
    controls: ol.control.defaults({
          attributionOptions: /** @type {olx.control.AttributionOptions} */ ({
               collapsible: false
          })
        }).extend([mousePositionControl]),
        layers: [
  new ol.layer.Tile({
            source: new ol.source.TileWMS({
      // MapServer URL this must be the same as the one declared in the
                      // MapServer .map file.
                          url : 'http://mapserver:8080/cgi-bin/mapserv.exe?map=wales_imagery.map',
  params : {
                                 // This must be the layer name as defined in the capabilities document
         'LAYERS' : 'wales_imagery',  
  'TILED': true,
'FORMAT': 'image/png',
'TRANSPARENT': 'true',
serverType: 'mapserver',
crossOrigin: 'anonymous'
  },
  projection : 'EPSG:4326',
  wrapX : false
}),
//opacity: 0.5
          })
        ],
        target: 'map',
        view: new ol.View({
          center: ol.proj.fromLonLat([-4.0, 51.7]),   // initial view position
          zoom: 9  // larger more zoomed in
        })
      });
  
      var projectionSelect = document.getElementById('projection');
      projectionSelect.addEventListener('change', function(event) {
        mousePositionControl.setProjection(event.target.value);
      });

      var precisionInput = document.getElementById('precision');
      precisionInput.addEventListener('change', function(event) {
        var format = ol.coordinate.createStringXY(event.target.valueAsNumber);
        mousePositionControl.setCoordinateFormat(format);
      });  
    </script>
  </body>
</html>


The bits in bold are what need to be configured correctly.

You can find out what layers are supported by a WMS Server using the GetCapabilities command:

http://127.0.0.1:8080/cgi-bin/mapserv.exe?map=wales_imagery.map&SERVICE=WMS&REQUEST=GetCapabilities

Make sure the ip address and port are correct.

What you are looking for is:

<Layer>
<Name>wales_imagery</Name>
<Title>WMS Server</Title>
<Abstract>This server was setup by damian.dixon at gmail.com</Abstract>
<CRS>EPSG:4326</CRS>

I was the Team Leader for MapLink Pro. Neil Tippett wrote a large part of the MapLink WMS Server. I ported the WMS server to Solaris and Linux and implemented the Xlib rendering and fine tuned it to run on multiple threads.

Friday, 15 September 2017

Debug output on Windows

On Windows you don't usually run anything from a console window. Even if you did printf/cout will not normally appear unless the application is a console application.

The following code can be used to out put to the Visual Studio debug console window.

It implements a variable argument list similar to printf.

#include <cstdarg>
#include <mutex>

std::mutex outputLock_;

void DebugMessage(const char * format, ...)
{
    char buffer[1024];
    va_list marker;
    va_start(markerformat);
    vsnprintf(buffersizeof(buffer), formatmarker);
    buffer[sizeof(buffer) - 1] = '\0';
 
    std::lock_guard<std::mutex> guard(outputLock_);
    OutputDebugStringA(buffer);
    va_end(marker);
}

A number of points to note:
  • It is thread safe but the output may be scrambled.
  • Limited buffer size, though this is safe.
  • Only supports code page output.
  • If you replace OutputDebugStringA with printf/cout/cerr the method becomes more portable.

Saturday, 2 September 2017

Orthorectification

Orthorectification is required to align imagery with a Map.

Normally this is done by the supplier of the imagery.

However imagery can be supplied without this process having been done.

In this case you have to do this yourself. There are quite a few Open Source tools that can perform this.

However you need the RCP data to do this or you need to create this information from the camera model.

This post is essentially to mark in time my thoughts so that I can retain IPR.

For UAV systems with cameras orthorectiification needs to occur in realtime to allow camera video or stills to overlay the mapping data correctly.

The simplest approaches are:

  • OpenGL using a "shadow mapping"/"image projection" approach to project the texture of each  frame onto a WGS84 ellipsoid with Terrain mesh. View looking straight down.
  • Using OpenCL/OpenGL to implement a reproection of the texture image.
While the above are not even 99% accurate the aim is to be close enough to achieve realtime display of the imagery on the map with a close enough image registration to the map.





Sunday, 20 August 2017

Converter programs for: CADRG, CIB, ECRG, GeoTIFF, RDTED

I've written software in the past to create CADRG/CIB and RDTED when I worked at Envitia on MapLink Pro. After I left Envitia I wrote a GeoTIFF exporter for a customer that needed to deploy GIS data to a MapServer instance.

CADRG/CIB is not particularly straight forward. The compression is non-standard and highly lossey. This is designed for display of 8bit imagery in a fast way.

ECRG is just a variation on CADRG but with JPEG2000 compression and an XML Table of Contents file. ECRG is a natural progression from CADRG/CIB and is for 24bit imagery.

Both ECRG and CADRG/CIB are Raster Product Formats (RPF) products which are wrapped in NITF.

GeoTIFF in comparison is straightforward.

If you require a tool that converts from common GIS formats to one of the following please contact me (damian.dixon at gmail.com) to discuss costs:

  • CADRG/CIB
  • ECRG
  • RDTED
  • GeoTIFF
The easiest support is for formats that are supported out of the box. For other formats there will be a Non-Recurring engineering cost.

For CADRG/CIB/ECRG/RDTED there will be at least 2 to 3 month lead time depending upon my availability. Once one of these formats is supported the lead time for each format will drop.

Glossary


  • CADRG - Compressed ARC Digitised Raster Graphics.
  • CIB - Controlled Image Base
  • ECRG - Enhanced Compressed ARC Raster 
  • RPF - Raster Product Format
  • NITF  - National Imagery Transfer Format
  • RDTED - DTED stored in NITF/RPF format.

Controlled Image Base and CIB are registered trademarks of the National Geospatial-Intelligence Agency, United States Department of Defense

Saturday, 19 August 2017

Reading Windows Bitmap data into an RGBA buffer.

To be able to save the image data into a TIFF tile we need the data in RGBA format. Windows provides the data in BGRA.

The following code deals with this:

    // copy result into buffer
    unsigned char *bits = Bits;
    for (uint32_t row = 0; row < (uint32_t)info.dsBm.bmHeight; ++row)
    {
        // This is windows the image is upside down.
        size_t line = (info.dsBm.bmHeight - 1) - row;
        unsigned char *pcopy_to = buffer + (line * stride);

        // Only copy stride as it may be smaller than scanline in bitmap.
        // Windows is also BGR rather than RGB
        uint8_t *pbits = bits;
        uint32_t *ptr = (uint32_t*)pcopy_to;
        uint32_t *ptrEnd = (uint32_t*)(pcopy_to + stride);
        while (ptr != ptrEnd)
        {
            *ptr = (0xff << 24) | (pbits[0] << 16) |
                     (pbits[1] << 8) | pbits[2];
            ++ptr;
            pbits += 4;
        }

        bits += info.dsBm.bmWidthBytes; // need to advance by bitmap width not buffer width
    }

Not the most elegant of code but its fast enough for what I currently need.

An extension to the above is to look for the RGB unique value set into the buffer before drawing (see previous post) ignoring the alpha value. This gets around issues when drawing rasters with masks using GDI.


Quick Clear of Windows 32bit Bitmap to RGBA

Clearing a 32bit Bitmap on Windows to a specific BGRA in a fairly quick way.

Get information about the Bitmap so that we can gain access to the image bits:

    HBITMAP hbmpOld = (HBITMAP)SelectObject(CompatibleHDC, Bitmap);
    DIBSECTION info = { 0 };
    ::GetObject(Bitmap, sizeof(DIBSECTION), &info);

Clear the Bitmap to a known BGRA:

    uint32_t value = (((uint8_t)255) << 24) |
                          (Parameters.BackgroundBlue << 16) | 
                          (Parameters.BackgroundGreen << 8) | 
                          (Parameters.BackgroundRed);
    unsigned char *bits = Bits;
    for (uint32_t row = 0; row < (uint32_t)info.dsBm.bmHeight; ++row)
    {
       uint32_t *to_set = (uint32_t *)(((uint8_t*)info.dsBm.bmBits) +
                               (row * info.dsBm.bmWidthBytes));
       uint32_t *ptrEnd = (uint32_t*)((uint8_t *)to_set +
                                            info.dsBm.bmWidthBytes);

       while (to_set != ptrEnd)
       {
         *to_set = value;
         ++to_set;
       }
    }

Alpha is being set to 255. All data is 4 byte aligned by definition on Windows in this case. Width of the bitmap is padded to 4 bytes and we've created a 32bit Bitmap.

Windows GDI only deals with BGR, unless you are using GDI+. GDI+ is however still not the answer for drawing with alpha on Windows, neither is Direct2D. A GDI draw will always set the alpha channel to 0.

The real difficult bit is picking an RGB value that does not appear in the output image you are drawing especially if you can't control the drawing.

You are going to have to post process the alpha channel after drawing on a per pixel basis.

MapServer and GeoTIFF mosaic example

I have been writing a GeoTIFF exporter from a propitiatory GIS system for a client who wishes to use MapServer as a WMS Server.

As part of the exporter I needed to write out an example MAP file which is pretty much complete with only a couple of modifications required by the user.

The following is an example of the MAP FILE:

MAP
 NAME "wales_imagery"   # This needs to be unique
 STATUS ON
 SIZE 400 300  # default size - kind of ignored for WMS
 # Symbol list - This is relative to the .map file.
 SYMBOLSET "../../apps/etc/symbols.txt"  # UPDATE
 IMAGECOLOR 255 255 255
 # Font list - This is relative to the.map file.
 FONTSET "../../apps/etc/fonts.txt"  # UPDATE

OUTPUTFORMAT
 NAME "png"
 DRIVER AGG/PNG
 MIMETYPE "image/png"
 IMAGEMODE RGBA
 TRANSPARENT ON
 EXTENSION "png"
 FORMATOPTION "GAMMA=0.75"
END
OUTPUTFORMAT
 NAME "gif"
 DRIVER GD/GIF
 MIMETYPE "image/gif"
 IMAGEMODE PC256
 EXTENSION "gif"
END
OUTPUTFORMAT
 NAME "png8"
 DRIVER AGG/PNG8
 MIMETYPE "image/png; mode=8bit"
 IMAGEMODE RGB
 EXTENSION "png"
 FORMATOPTION "QUANTIZE_FORCE=on"
 FORMATOPTION "QUANTIZE_COLORS=256"
 FORMATOPTION "GAMMA=0.75"
END
OUTPUTFORMAT
 NAME "jpeg"
 DRIVER AGG/JPEG
 MIMETYPE "image/jpeg"
 IMAGEMODE RGBA
 TRANSPARENT ON
 EXTENSION "jpg"
 FORMATOPTION "GAMMA=0.75"
END
OUTPUTFORMAT
 NAME "svg"
 DRIVER CAIRO/SVG
 MIMETYPE "image/svg+xml"
 IMAGEMODE RGBA
 EXTENSION "svg"
END
OUTPUTFORMAT
 NAME "pdf"
 DRIVER CAIRO/PDF
 MIMETYPE "application/x-pdf"
 IMAGEMODE RGBA
 TRANSPARENT ON
 EXTENSION "pdf"
END
OUTPUTFORMAT
 NAME "GTiff"
 DRIVER GDAL/GTiff
 MIMETYPE "image/tiff"
 IMAGEMODE RGB
 EXTENSION "tif"
END
OUTPUTFORMAT
 NAME "kml"
 DRIVER KML
 MIMETYPE "application/vnd.google-earth.kml.xml"
 IMAGEMODE RGBA
 TRANSPARENT ON
 EXTENSION "kml"
END
OUTPUTFORMAT
 NAME "kmz"
 DRIVER KMZ
 MIMETYPE "application/vnd.google-earth.kmz"
 IMAGEMODE RGBA
 TRANSPARENT ON
 EXTENSION "kmz"
END
OUTPUTFORMAT
 NAME "cairopng"
 DRIVER CAIRO/PNG
 MIMETYPE "image/png"
 IMAGEMODE RGBA
 TRANSPARENT ON
 EXTENSION "png"
END
OUTPUTFORMAT
 NAME "myUTFGrid"
 DRIVER UTFGRID
 FORMATOPTION "LABELS=true"
 FORMATOPTION "UTFRESOLUTION=4"
 FORMATOPTION "DUPLICATES=false"
END
#
# End of OUTPUTFORMAT definitions
#

 # Extents :
 EXTENT 304430.5800000000 5586579.0200000005 500356.8900000000 5874520.2599999998
 UNITS METERS # units the extent is in.

 # Where to find shapefile data - This is relative to the MapServer '.map' file
 SHAPEPATH "../../apps/maps/wales_imagery" # UPDATE

 # Maximum image request size
 MAXSIZE 4096

 WEB
# Where the GeoTIFF data can be found.This is relative to the MapServer '.map' file.
IMAGEPATH "../../apps/maps/wales_imagery" # UPDATE
IMAGEURL "/ms_tmp/"
METADATA
 "wms_title" "WMS Server"
 "wms_onlineresource" "http://mapserver:8080/cgi-bin/mapserv.exe?map=wales_imagery.map&" # UPDATE IP address/name

 # The projections offered via the WMS server
 "wms_srs" "EPSG:4326 EPSG:32630"  # UPDATE with any additional EPSG codes you want to offer.

 "wms_feature_info_mime_type" "text/plain"
 "wms_abstract" "This server was setup by damian.dixon at gmail.com"
 "ows_enable_request" "*"
END # Metadata
 END # WEB

 # Default EPSG of the data - MapServer will re-project the data to one of the requested 'wms_srs' EPSGs.
 PROJECTION
"init=EPSG:32630"
 END

 # Start of layer definitions
 #
 # Layers are defined in order of display(if not via a WMS)
 #
 LAYER
# WMS Layer name
NAME "wales_imagery"

# Metadata to return in capabilities document.
METADATA
 "wms_title" "wales_imagery"
 "wms_include_items" "all"
END

# Should be ON / DEFAULT / OFF
STATUS ON

# The projection the data is in.
PROJECTION
 # This is the data EPSG code.
 "init=EPSG:32630"
END

  # Specify the RGBA bands to use
  PROCESSING "BANDS=1,2,3,4"

# Specify the Shapefile that contains the tile index data.
TILEINDEX "wales_imagery.shp"
TILEITEM "Location"

# What the layer type is.
TYPE RASTER
 END # layer

END # Map File

GeoTIFF setup

The GeoTIFF is written out as a pyramid file with an Alpha band.

Because  the extent of the data may be quite huge the output is split into an 4 by 2 grid. Basically 8 top level GeoTIFF files.

The number of levels contained in the GeoTIFF is based on the data in the original GIS system and the number of levels of detail being exported.

If the map scale gap between the levels of detail is too great another pyramid level is added.

Each top level tile has a slight overlap to ensure that the images appear to be seamless.

Transparency

There are a couple of interesting points in the MAP file:
  1. Enabling Transparency in the WMS tile returned.
  2. Specifying a Transparent pixel in the input GeoTIFF.

Transparency in a WMS GetMap response

This is enabled by specifying the format is ARGB

OUTPUTFORMAT
  NAME "png"
  DRIVER AGG/PNG
  MIMETYPE "image/png"
  IMAGEMODE RGBA
  TRANSPARENT ON
  EXTENSION "png"
  FORMATOPTION "GAMMA=0.75"
END

Transparent pixel

There are two ways of specifying Alpha value.

The first is to specified in the LAYER section an RGB to convert to transparent as follows:

        OFFSET 255 255 255

The values specified represent the pixel colour value. White is a poor choice. Black can also be a poor choice as lettering is usually black. The colour needs to really be unique, i.e. not used in the image other than to indicate to mask the value out.

The second approach is to specify the alpha channel in the data as follows:

PROCESSING "BANDS=1,2,3,4"

The channel numbers represent R, G, B and alpha in that order.

I am using the second approach because I can generate an alpha channel. However to generate the alpha channel I am using exactly the approach that MapServer is using in the first approach. While this takes time to process you do save time and disk space when compressing and writing out the tile.

The following links have further information about transparency:

Combining the GeoTIFF Mosaic

The GeoTIFFs are combined using a Tiled Index using the following command:

gdaltindex wales_imagery.shp *.tif

The paths specified to the GeoTIFF files must be relative to the shapefile that is created and to where the GIS data will be deployed in MapServer.

In the MAP file I specify where the Shapefile and GeoTIFFs can be found:

SHAPEPATH "../../apps/maps/wales_imagery"
IMAGEPATH "../../apps/maps/wales_imagery"

In the LAYER section you specify you are using a Tile Index as follows:

# Specify the Shapefile that contains the tile index data.
TILEINDEX "wales_imagery.shp"
TILEITEM "Location"

You specify the shapefile as the TILEINDEX.

Deploying the MAP File

The quickest and easiest way to deploy the MAP file is to put it into the cgi-bin directory.

This will be immediately picked up by the Server and accessible. You can tweak the file at anytime and the changes are automatically picked up.

The data is normally specified relative to the MAP file.

Please contact me if you need a similar tool at damian.dixon at gmail.com to discuss.