Overlaying Multiple Layers

Creating a Composite Map of Multiple Layers from One Server

The WMS specification allows the server to merge multiple layers into a single raster map. Metacarta's 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')

Combining Layers from One Server with Data from Other Sources

This example, a continuation of the one preceding it, illustrates how you can merge the boundaries raster map with vector data. Combine two separate calls to vec2mtx to create a 4-color raster map showing interior land areas, coastlines, national boundaries, oceans, and world rivers. The vec2mtx function converts latitude-longitude vectors to a regular data grid.

  1. Read the landareas vector shapefile and convert it to an image.

    land = shaperead('landareas', 'UseGeoCoords', true);
    lat = [land.Lat];
    lon = [land.Lon];
    [landImage, refvec] = ...
       vec2mtx(lat, lon, 2, [-90, 90], [-180,180], 'filled');
    mergedImage = landImage;
  2. Read the worldrivers vector shapefile and convert it to an image.

    rivers = shaperead('worldrivers.shp','UseGeoCoords',true);
    riverImage = vec2mtx([rivers.Lat], [rivers.Lon], landImage, refvec);
  3. Merge the rivers with the land.

    mergedImage(riverImage == 1) = 3;
  4. Obtain the boundaries image from the VMAP0 server.

    vmap0 = wmsfind('vmap0.tiles', 'SearchField', 'serverurl');
    vmap0 = wmsupdate(vmap0);
    layer = refine(vmap0, 'country_01');
    height = size(riverImage, 1);
    width = size(riverImage, 2);
    [boundaries, R] = wmsread(layer, 'ImageFormat', 'image/png', ...
        'ImageHeight', height, 'ImageWidth', width);
  5. Merge the rivers and land with the boundaries. You must first flip mergedImage because it has a reference defined by a referencing vector, where columns run from south to north. (The columns of WMS images run from north to south.)

    mergedImage = flipud(mergedImage);
    index = boundaries(:,:,1) ~= 255 ...
        & boundaries(:,:,2) ~= 255 ...
        & boundaries(:,:,3) ~= 255;
    mergedImage(index) = 1;
  6. Display the result.

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

Drape Topography and Ortho-Imagery Layers over a Digital Elevation Model Layer

  1. Uncompress the zip file and read it with the usgs24kdem function. Set the geographic limits to the minimum and maximum values in the DEM file.

    filenames = gunzip('sanfranciscos.dem.gz', tempdir);
    demFilename = filenames{1};
    [lat,lon,Z,header,profile] = usgs24kdem(demFilename, 1);
    Z(Z==0) = -1;
    latlim = [min(lat(:)) max(lat(:))];
    lonlim = [min(lon(:)) max(lon(:))];
  2. Display the USGS 24K DEM data. Create map axes for the United States and assign an appropriate elevation colormap. Use daspectm to display the elevation data.

    usamap(latlim, lonlim)
    geoshow(lat, lon, Z, 'DisplayType','surface')
    title('San Francisco South 24K DEM');

  3. Set the point of view by adjusting the CameraPosition, CameraTarget, and CameraAngle axes properties.

    cameraPosition = [87907 4.4313e+06 53103];
    cameraTarget = [-1961.5 4.492e+06 231.46];
    cameraAngle = 5.57428;
    set(gca,'CameraPosition', cameraPosition, ...
        'CameraTarget', cameraTarget, ...
        'CameraViewAngle', cameraAngle)

  4. The USGS National Map provides ortho-imagery and topographic maps from various regions of the United States. One method to obtain the high-resolution ortho-imagery layer name is to obtain the capabilities document from the server. The ortho-imagery layer is the only layer from this 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)
  5. Use the USGS topographic base map layer.

    topoLayer = wmsfind('usgs*topographic*small*scale','SearchField','any');
  6. Set the size for the requested images.

    imageHeight = size(Z,1)*3;
    imageWidth  = size(Z,2)*3;
  7. Request a map of the ortho-imagery layer.

    orthoImage = wmsread(orthoLayer, 'Latlim', latlim, 'Lonlim', lonlim, ...
        'ImageHeight', imageHeight, 'ImageWidth',  imageWidth);
  8. Request a map of the topographic layer.

    topoMap = wmsread(topoLayer, 'Latlim', latlim, 'Lonlim', lonlim, ...
        'ImageHeight', imageHeight, 'ImageWidth',  imageWidth);
  9. Drape the ortho-image onto the elevation data.

    usamap(latlim, lonlim)
    geoshow(lat, lon, Z, ...
        'DisplayType', 'surface', 'CData', orthoImage);
    title('San Francisco Ortho-Image');
    axis vis3d
    set(gca,'CameraPosition', cameraPosition, ...
        'CameraTarget', cameraTarget, ...
        'CameraViewAngle', cameraAngle)

  10. Drape the topographic map onto the elevation data.

    usamap(latlim, lonlim)
    geoshow(lat, lon, Z, ...
        'DisplayType', 'surface', 'CData', topoMap);
    title('San Francisco Topo Map');
    axis vis3d
    set(gca,'CameraPosition', cameraPosition, ...
        'CameraTarget', cameraTarget, ...
        'CameraViewAngle', cameraAngle)

Was this topic helpful?