Read WMS Maps Using Different Coordinate Reference Systems

To read a WMS map in the EPSG:4326 coordinate reference system, use the wmsread function. EPSG:4326 is based on the 1984 World Geodetic System (WGS84) datum. All servers in the WMS Database, and presumably all WMS servers in general, use the EPSG:4326 coordinate reference system. This system is a requirement of the OGC® WMS specification.

To read a WMS map in a different coordinate reference system, create a request URL using WebMapServer and WMSMapRequest objects and then read the map using the getMap method. For more information about coordinate reference systems, see Spatial Reference.

This example shows how to read a WMS map with data in Web Mercator coordinates, also known as WGS 84/Pseudo-Mercator coordinates. The Web Mercator coordinate system is commonly used by web applications.

Find Web Mercator Coordinates for Region

Import a GeoTIFF image of Boston as an array and a MapCellsReference object. Find the projected coordinate reference system for the image.

[A,R] = readgeoraster('boston.tif');
p_geotiff = R.ProjectedCRS;

Unproject the x- and y-world limits and then reproject them to use the Web Mercator coordinate reference system (EPSG:3857).

[latlim,lonlim] = projinv(p_geotiff,R.XWorldLimits,R.YWorldLimits);
p_webmercator = projcrs(3857);
[xlimits,ylimits] = projfwd(p_webmercator,latlim,lonlim);

To obtain imagery in this coordinate reference system, you must use WMSMapRequest and WebMapServer objects because the wmsread function only reads data in the WGS84 coordinate reference system (EPSG:4326).

Read and Display Map for Region

The USGS National Map provides orthoimagery and topography maps for various regions of the United States. The USGS Imagery Only server provides data in both WGS84 coordinates and Web Mercator coordinates.

Search the WMS Database for the USGS Imagery Only server and select the first layer.

doqLayer = wmsfind('usgsimageryonly','SearchField','serverurl');
doqLayer = wmsupdate(doqLayer);

Create WebMapServer and WMSMapRequest objects.

server = WebMapServer(doqLayer.ServerURL);
request = WMSMapRequest(doqLayer,server);

Modify the map request by setting properties of the WMSMapRequest object. For this example, specify an image height and width for a sample size of 5 meters. Set the map limits to cover the same region as the GeoTIFF file.

metersPerSample = 5;
h = round(diff(ylimits)/metersPerSample);
w = round(diff(xlimits)/metersPerSample);

ylimits = [ylimits(1), ylimits(1) + h*metersPerSample];
xlimits = [xlimits(1), xlimits(1) + w*metersPerSample];

request.CoordRefSysCode = 'EPSG:3857';
request.ImageHeight = h;
request.ImageWidth  = w;
request.XLim = xlimits;
request.YLim = ylimits;

Read a map of the orthoimagery in Web Mercator coordinates.

A_webmercator = getMap(server,request.RequestURL);
R_webmercator = request.RasterReference;

Display the orthoimagery on a map.

figure
mapshow(A_webmercator,R_webmercator)
axis tight
title({'USGS Digital Ortho-Quadrangle - Boston','Web Mercator'})

Read Boston place names from a shapefile. Unproject and reproject the place names so that they are in Web Mercator coordinates. Display the place names on the map.

S = shaperead('boston_placenames.shp');
x_names = [S.X] * unitsratio('sf','meter');
y_names = [S.Y] * unitsratio('sf','meter');
names = {S.NAME};

[lat_placenames,lon_placenames] = projinv(p_geotiff, ...
                                          x_names,y_names);
[x_webmercator,y_webmercator] = projfwd(p_webmercator, ...
                                        lat_placenames,lon_placenames);

text(x_webmercator,y_webmercator,names, ...
    'BackgroundColor',[0.9 0.9 0],'FontSize',6,'Clipping','on')

Compare the map you read from the server to the GeoTIFF image.

figure
mapshow('boston.tif')
axis tight
title({'boston.tif', p_geotiff.Name})
text(x_names,y_names,names, ...
    'BackgroundColor',[0.9 0.9 0],'FontSize',6,'Clipping','on')

You can also compare the maps by using geographic axes and a satellite basemap. When you plot on geographic axes, coordinates are referenced to the WGS84 coordinate reference system.

figure
geolimits(latlim,lonlim)
geobasemap('satellite')
text(lat_placenames,lon_placenames,names, ...
    'BackgroundColor',[0.9 0.9 0],'FontSize',6)
title({'Satellite Basemap','WGS84'})