Наложение нескольких слоев

Создайте составную карту нескольких слоев из одного сервера

Спецификация WMS позволяет серверу объединять несколько слоев в одну растровую карту. Сервер VMAP0 MetaCarta содержит много слоев данных, таких как береговые линии, национальные контуры, океан и земля. Чтение и отображение композита из нескольких слоев с сервера VMAP0. Визуализированная карта имеет пространственное разрешение 0,5 степень на камеру.

  1. Найдите и обновите слои VMAP0.

    vmap0 = wmsfind('vmap0.tiles', 'SearchField', 'serverurl');
    vmap0 = wmsupdate(vmap0);
    
  2. Создайте массив из нескольких слоев, которые включают землю и океан, береговые линии и национальные контуры.

    layers = [refine(vmap0, 'coastline_01'); ...
              refine(vmap0, 'country_01'); ...
              refine(vmap0, 'ground_01'); ...
              refine(vmap0, 'inwater'); ...
              refine(vmap0, 'ocean')];
    
  3. Найдите составную карту. Запросите размер камеры 5 степеней путем установки параметров высоты изображения и ширины изображения. Установите значение 'Transparent' на true, чтобы все пиксели, не представляющие функции или значения данных в слое, были установлены на прозрачное значение в полученном изображении, что позволяет создать составную карту.

    [overlayImage, R] = wmsread(layers, 'Transparent', true, ...
                                'ImageHeight', 360, 'ImageWidth', 720);
    
  4. Отобразите композитную карту.

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

Объединение слоев с одного сервера с данными из других источников

Этот пример иллюстрирует, как можно объединить растровую карту контуров с векторными данными.

  1. Размещение глобального растра с камерами 1/2 градуса с помощью georefcells функция. Задайте столбцы с севера на юг для согласованности с wmsread. Результат R является географическим объектом растровой привязки.

    latlim = [-90 90];
    lonlim = [-180 180];
    cellExtent = 1/2;
    R = georefcells(latlim,lonlim,  ...
        cellExtent,cellExtent,'ColumnsStartFrom','north');
    
  2. Чтение landareas polygon shapefile и преобразовать его в растровую карту.

    land = shaperead('landareas.shp', 'UseGeoCoords', true);
    lat = [land.Lat];
    lon = [land.Lon];
    land = vec2mtx(lat, lon, zeros(R.RasterSize),R, 'filled');
    
  3. Чтение worldrivers polyline shapefile и преобразуйте его в растровую карту.

    riverLines = shaperead('worldrivers.shp','UseGeoCoords',true);
    rivers = vec2mtx([riverLines.Lat], [riverLines.Lon], land, R);
    
  4. Слить реки с землей.

    merged = land;
    merged(rivers == 1) = 3;
    
  5. Получите информацию о системе координат-привязок из landareas shapefile с использованием shapeinfo функция. The worldrivers shapefile использует ту же систему координат-привязок. Установите GeographicCRS свойство объекта ссылки.

    info = shapeinfo('landareas.shp');
    R.GeographicCRS = info.CoordinateReferenceSystem;
  6. Получите изображение контуров от сервера VMAP0.

    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);
    
  7. Подтвердите, что контуры и объединенные растры совпадают.

    isequal(boundariesR,R)
    
    ans =
    
      logical
    
       1
  8. Объедините реки и сушу с контурами.

    index = boundaries(:,:,1) ~= 255 ...
        & boundaries(:,:,2) ~= 255 ...
        & boundaries(:,:,3) ~= 255;
    merged(index) = 1;
    
  9. Отобразите результат.

    figure
    worldmap(merged, R)
    geoshow(merged, R, 'DisplayType', 'texturemap')
    colormap([.45 .60 .30; 0 0 0; 0 0.5 1; 0 0 1])

Drape ортоимагерия над DEM

Считайте данные о повышении и ссылку на географические проводки для области вокруг Саут-Боулдер-Пик в Колорадо. Затем обрезайте данные повышения на меньшую площадь, используя 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);

Отображение данных о повышении. Для этого создайте набор осей карты для США, постройте график данных как поверхности и примените соответствующую палитру. Просмотрите карту в 3-D путем настройки положения камеры и цели. Установите вертикальное преувеличение с помощью 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)

Перетащите ортоизображение над данными о повышении. Для этого сначала получите имена слоев ортоимагерии высокого разрешения из Национальной карты USGS с помощью wmsinfo. В этом случае слой ортоимагерии является единственным слоем от сервера. Используйте несколько попыток подключиться к серверу, если он занят.

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

Запросите карту слоя ортоимагерии, используя wmsread. Чтобы отобразить ортоимагерию, используйте geoshow и установите CData свойство слою.

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)

Файл DTED, используемый в этом примере, предоставлен Геологической службой США.

См. также

| |

Похожие темы