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. 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')

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 out a global raster with 1/2-degree cells. Specify columns running north-to-south, for consistency with wmsread.

    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', '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. 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);
  6. Confirm that the boundaries and merged rasters are coincident.

  7. Merge the rivers and land with the boundaries.

    index = boundaries(:,:,1) ~= 255 ...
        & boundaries(:,:,2) ~= 255 ...
        & boundaries(:,:,3) ~= 255;
    merged(index) = 1;
  8. 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 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('/USGSImageryTopo/','SearchField','serverurl');
  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)

See Also

| |

Related Topics