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.

Get the capabilities document for the Metacarta VMAP0 server. Get the list of layers from the capabilities document.

info = wmsinfo('http://vmap0.tiles.osgeo.org/wms/vmap0');
vmap0 = info.Layer;

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 0.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 the composite map.

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

The data used in this example is from Metacarta.

Combine Layers from One Server with Data from Other Sources

This example shows how to merge a boundaries raster map with vector data.

Specify a global raster with 0.5-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 = 0.5;
R = georefcells(latlim,lonlim,  ...
    cellExtent,cellExtent,'ColumnsStartFrom','north');

Read a shapefile that contains global land area polygons 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');

Read a shapefile that contains world river polylines 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;

Get coordinate reference system information from the land areas shapefile by using the shapeinfo function. The world rivers shapefile uses the same coordinate reference system. Set the GeographicCRS property of the reference object.

info = shapeinfo('landareas.shp');
R.GeographicCRS = info.CoordinateReferenceSystem;

Read the boundaries image from the VMAP0 server.

info = wmsinfo('http://vmap0.tiles.osgeo.org/wms/vmap0');
vmap0 = info.Layer;
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)
ans = logical
   1

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([0.45 0.60 0.30; 0 0 0; 0 0.5 1; 0 0 1])

The data used in this example is from the US National Geospatial-Intelligence Agency (NGA) and Metacarta.

Drape Orthoimagery Over DEM

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

[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 an axesm-based map 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 by using the daspectm function.

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 the wmsinfo function. 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 the wmsread function. To display the orthoimagery, use the geoshow function 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 from the US Geological Survey.

See Also

| |

Related Topics