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.

Find and update the VMAP0 layers.

vmap0 = wmsfind('vmap0.tiles', 'SearchField', 'serverurl');
vmap0 = wmsupdate(vmap0);

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')];

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);
```

Display your composite map.

figure
worldmap('world')
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.

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,  ...
    cellExtent,cellExtent,'ColumnsStartFrom','north')

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

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);

Merge the rivers with the land.

merged = land;
merged(rivers == 1) = 3;

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);

Confirm that the boundaries and merged rasters are coincident.
```
isequal(boundariesR,R)
```

Merge the rivers and land with the boundaries.

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

Display the result.

figure
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

Open Live Script

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.

figure
usamap(R.LatitudeLimits, R.LongitudeLimits)
geoshow(Z,R,'DisplayType','surface')
demcmap(Z)
title('Elevation');

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

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 = ...
   'http://basemap.nationalmap.gov/ArcGIS/services/USGSImageryOnly/MapServer/WMSServer?';
while(isempty(info))
    try
        info = wmsinfo(serverURL);
        orthoLayer = info.Layer(1);
    catch e         
        attempt = attempt + 1;
        if attempt > numberOfAttempts
            throw(e);
        else
            fprintf('Attempting to connect to server:\n"%s"\n', serverURL)
        end        
    end
end

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);

figure
usamap(R.LatitudeLimits,R.LongitudeLimits)
geoshow(Z,R,'DisplayType','surface','CData',orthoImage);
title('Orthoimage Draped Over Elevation');
set(gca,'CameraPosition', cameraPosition, ...
        'CameraTarget', cameraTarget)
daspectm('m',3)

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

Documentation