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.