Main Content

Overlay Multiple Layers

Create Composite Map of Multiple Layers from One Server

The WMS specification allows the server to merge multiple layers into a single raster map. The MetaCarta VMAP0 server contains many data layers, such as coastlines, national boundaries, ocean, and ground. Read and display a composite of multiple layers from the VMAP0 server. The rendered map has a spatial resolution of 0.5 degrees per cell.

  1. Find and update the VMAP0 layers.

    vmap0 = wmsfind('vmap0.tiles', 'SearchField', 'serverurl');
    vmap0 = wmsupdate(vmap0);
  2. Create an array of multiple layers that include ground and ocean, coastlines and national boundaries.

    layers = [refine(vmap0, 'coastline_01'); ...
              refine(vmap0, 'country_01'); ...
              refine(vmap0, 'ground_01'); ...
              refine(vmap0, 'inwater'); ...
              refine(vmap0, 'ocean')];
  3. Retrieve the composite map. Request a cell size of .5 degrees by setting the image height and image width parameters. Set 'Transparent' to true so that all pixels not representing features or data values in a layer are set to a transparent value in the resulting image, making it possible to produce a composite map.

    [overlayImage, R] = wmsread(layers, 'Transparent', true, ...
                                'ImageHeight', 360, 'ImageWidth', 720);
  4. Display your composite map.

    geoshow(overlayImage, R);
    title('Composite of VMAP0 Layers')

Combine Layers from One Server with Data from Other Sources

This example illustrates how you can merge a boundaries raster map with vector data.

  1. Layout a global raster with 1/2-degree cells by using the georefcells function. Specify columns running north-to-south, for consistency with wmsread. The result R is a geographic raster reference object.

    latlim = [-90 90];
    lonlim = [-180 180];
    cellExtent = 1/2;
    R = georefcells(latlim,lonlim,  ...
  2. Read the landareas polygon shapefile and convert it to a raster map.

    land = shaperead('landareas.shp', 'UseGeoCoords', true);
    lat = [land.Lat];
    lon = [land.Lon];
    land = vec2mtx(lat, lon, zeros(R.RasterSize),R, 'filled');
  3. Read the worldrivers polyline shapefile and convert it to a raster map.

    riverLines = shaperead('worldrivers.shp','UseGeoCoords',true);
    rivers = vec2mtx([riverLines.Lat], [riverLines.Lon], land, R);
  4. Merge the rivers with the land.

    merged = land;
    merged(rivers == 1) = 3;
  5. Get coordinate reference system information from the landareas shapefile using the shapeinfo function. The worldrivers shapefile uses the same coordinate reference system. Set the GeographicCRS property of the reference object.

    info = shapeinfo('landareas.shp');
    R.GeographicCRS = info.CoordinateReferenceSystem;
  6. Obtain the boundaries image from the VMAP0 server.

    vmap0 = wmsfind('vmap0.tiles', 'SearchField', 'serverurl');
    vmap0 = wmsupdate(vmap0);
    layer = refine(vmap0, 'country_01');
    height = R.RasterSize(1);
    width  = R.RasterSize(2);
    [boundaries,boundariesR] = wmsread(layer, 'ImageFormat', 'image/png', ...
        'ImageHeight', height, 'ImageWidth', width);
  7. Confirm that the boundaries and merged rasters are coincident.

    ans =
  8. Merge the rivers and land with the boundaries.

    index = boundaries(:,:,1) ~= 255 ...
        & boundaries(:,:,2) ~= 255 ...
        & boundaries(:,:,3) ~= 255;
    merged(index) = 1;
  9. Display the result.

    worldmap(merged, R)
    geoshow(merged, R, 'DisplayType', 'texturemap')
    colormap([.45 .60 .30; 0 0 0; 0 0.5 1; 0 0 1])

Drape Orthoimagery Over DEM

Read elevation data and a geographic postings reference for an area around South Boulder Peak in Colorado. Then, crop the elevation data to a smaller area using geocrop.

[fullZ,fullR] = readgeoraster('n39_w106_3arc_v2.dt1','OutputType','double');

latlim = [39.25 40.0];
lonlim = [-106 -105.5];
[Z,R] = geocrop(fullZ,fullR,latlim,lonlim);

Display the elevation data. To do this, create a set of map axes for the United States, plot the data as a surface, and apply an appropriate colormap. View the map in 3-D by adjusting the camera position and target. Set the vertical exaggeration using daspectm.

usamap(R.LatitudeLimits, R.LongitudeLimits)

cameraPosition = [218100 4367600 183700];
cameraTarget = [0 4754200 2500];
set(gca,'CameraPosition', cameraPosition, ...
        'CameraTarget', cameraTarget)

Drape an orthoimage over the elevation data. To do this, first get the names of high-resolution orthoimagery layers from the USGS National Map using wmsinfo. In this case, the orthoimagery layer is the only layer from the server. Use multiple attempts to connect to the server in case it is busy.

numberOfAttempts = 5;
attempt = 0;
info = [];
serverURL = ...
        info = wmsinfo(serverURL);
        orthoLayer = info.Layer(1);
    catch e         
        attempt = attempt + 1;
        if attempt > numberOfAttempts
            fprintf('Attempting to connect to server:\n"%s"\n', serverURL)

Request a map of the orthoimagery layer using wmsread. To display the orthoimagery, use geoshow and set the CData property to the layer.

imageHeight = size(Z,1);
imageWidth  = size(Z,2);

orthoImage = wmsread(orthoLayer,'Latlim',R.LatitudeLimits, ...
    'Lonlim',R.LongitudeLimits,'ImageHeight', imageHeight, ...
    'ImageWidth',  imageWidth);

title('Orthoimage Draped Over Elevation');
set(gca,'CameraPosition', cameraPosition, ...
        'CameraTarget', cameraTarget)

The DTED file used in this example is courtesy of the US Geological Survey.

See Also

| |

Related Topics