does anyone know how to extract ONLY the elevation of a country ?

I am capable to use this

http://www.mathworks.com/help/map/retrieving-elevation-data.html

however , i would like to have ONLY the elevation of a country and not of the surrounding countries....

e.g. for cameroon ...

I tried to trim the elevation along the border of a country (imported as shapefile) but I am not really successful ...

please help

THX

M

[Matteo]

btw, the code I am using is the following

close all clear all

layers = wmsfind('nasa.network*elev', 'SearchField', 'serverurl'); layers = wmsupdate(layers);

disp(layers,'Properties',{'LayerTitle','LayerName'})

aster = layers.refine('earthaster', 'SearchField', 'layername');

latlim = [0 15]; lonlim = [5 20];

cellSize = dms2degrees([0,1,0]); [ZA, RA] = wmsread(aster, 'Latlim', latlim, 'Lonlim', lonlim, ... 'CellSize', cellSize, 'ImageFormat', 'image/bil');

figure('Renderer','zbuffer') worldmap('Cameroon') geoshow(ZA, RA, 'DisplayType', 'texturemap') demcmap(double(ZA)) title({'Cameroon and Surrounding Region', aster.LayerTitle});

vmap0 = wmsfind('vmap0.tiles', 'SearchField', 'serverurl'); boundaries = refine(vmap0, 'country_02'); B = wmsread(boundaries, 'Latlim', latlim, ... 'Lonlim', lonlim, 'CellSize', cellSize, 'ImageFormat','image/png'); ZB = ZA; ZB(B(:,:,1) < 250) = min(ZA(:));

figure('Renderer','zbuffer') worldmap('Cameroon') demcmap(double(ZA)) geoshow(ZB, RA, 'DisplayType', 'texturemap') title({'Cameroon and Country Boundaries', aster.LayerTitle});

%%%%%%%%

worldmap({'Cameroon'}) land = shaperead('landareas.shp', 'UseGeoCoords', true); geoshow(land, 'FaceColor', [0.15 0.5 0.15]) cities = shaperead('worldcities', 'UseGeoCoords', true); geoshow(cities, 'Marker', '.', 'Color', 'red')

%%%%%%%

figure worldmap 'cameroon' axis off getm(gca,'MapProjection')

geoshow('landareas.shp', 'FaceColor', [0.5 0.7 0.5]) geoshow('worldrivers.shp', 'Color', 'blue') geoshow('worldcities.shp', 'Marker', '.', 'Color', 'red')

%%%%%%%

%% load files % % url = 'http://www.zmuc.dk/public/GIS-intro/GIS-data/Lektion1/Temaer/'; % % try % % urlwrite([url '/cntry98.shp'],'cntry98.shp'); % % urlwrite([url '/cntry98.dbf'],'cntry98.dbf'); % % urlwrite([url '/cntry98.shx'],'cntry98.shx'); % % catch % % error('Unable to download cntry98 shapefiles.'); % % end % % % % % % % % % % % % url = 'http://openmap.bbn.com/data/shape/cntry02.tar.gz'; % % try % % files = untar(url) % % catch % % error('Unable to download %s.', url); % % end % % % % % % % % % % % % url = 'http://mappinghacks.com/data/world_borders.zip'; % % try % % files = unzip(url) % % catch % % error('Unable to download %s.', url); % % end

figure worlddatamap cntry98.shp world title({'World Countries','Using cntry98.shp'})

figure worlddatamap('cntry98.shp','cameroon','line')

figure worlddatamap('cntry98.shp','cameroononly','line')

figure worlddatamap('cntry02.shp','cameroon','cities') %title({'Cities of Cameroon','Using cntry02.shp'}) title({'Cities of Cameroon'})

%%%%%%%%%%% elevation and boundaries

layers = wmsfind('nasa.network*elev', 'SearchField', 'serverurl'); layers = wmsupdate(layers);

disp(layers,'Properties',{'LayerTitle','LayerName'})

aster = layers.refine('earthaster', 'SearchField', 'layername');

latlim = [0 15]; lonlim = [5 20];

cellSize = dms2degrees([0,1,0]); [ZA, RA] = wmsread(aster, 'Latlim', latlim, 'Lonlim', lonlim, ... 'CellSize', cellSize, 'ImageFormat', 'image/bil');

figure('Renderer','zbuffer') worldmap('Cameroon') geoshow(ZA, RA, 'DisplayType', 'texturemap') demcmap(double(ZA)) title({'Cameroon and Surrounding Region', aster.LayerTitle});

hold on worlddatamap('cntry02.shp','cameroon','cities') %title({'Cities of Cameroon','Using cntry02.shp'}) title({'Cities of Cameroon'})

%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%% PLOT CORRECT %%%%%

hold on worlddatamap('cntry02.shp','cameroon','cities') %title({'Cities of Cameroon','Using cntry02.shp'}) title({'Cities of Cameroon'})

figure('Renderer','zbuffer') worldmap('Cameroon') geoshow(ZA, RA, 'DisplayType', 'texturemap') demcmap(double(ZA)) title({'Cameroon and Surrounding Region', aster.LayerTitle});

vmap0 = wmsfind('vmap0.tiles', 'SearchField', 'serverurl'); boundaries = refine(vmap0, 'country_02'); B = wmsread(boundaries, 'Latlim', latlim, ... 'Lonlim', lonlim, 'CellSize', cellSize, 'ImageFormat','image/png'); ZB = ZA; ZB(B(:,:,1) < 250) = min(ZA(:));

figure('Renderer','zbuffer') worldmap('Cameroon') demcmap(double(ZA)) geoshow(ZB, RA, 'DisplayType', 'texturemap') title({'Cameroon and Country Boundaries', aster.LayerTitle});

figure('Renderer','zbuffer') worldmap('Cameroon') imagesc((ZA)) geoshow(ZB, RA, 'DisplayType', 'texturemap') title({'Cameroon and Country Boundaries', aster.LayerTitle});

%%%

figure,

imagesc(ZA)

%%%

index = ZA < 0;

ZA(index)=0;

dlmwrite('ZAzero.mat',ZA,'precision',6,'delimiter',' ')

% % ZA(ZA==0)= NaN ; % figure % imagescnan(ZA) %

% basename = 'ZAzero'; % imagefile = [basename '.jpg']; % RGB = imread(imagefile); % worldfile = getworldfilename(imagefile); % R = worldfileread(worldfile, 'geographic', size(RGB)); % filename = [basename '.tif']; % geotiffwrite(filename, RGB, R) % figure % usamap(RGB, R) % geoshow(filename) %

M.grid = ZA; M.dx = 1; M.dy = 1;

M

% matteo=[10 11 24 2 5 10; 5 10 14 21 5 13; 11 31 14 2 5 1 ] % % matteogr.grid = matteo; % matteogr.dx = 1; % matteogr.dy = 1; % % matteogr

%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%% PLOT CORRECT ONLY CAMEROON %%%%%

hold on worlddatamap('cntry02.shp','cameroononly','cities') %title({'Cities of Cameroon','Using cntry02.shp'}) title({'Cities of Cameroon'})

figure('Renderer','zbuffer') worldmap('Cameroon') geoshow(ZA, RA, 'DisplayType', 'texturemap') demcmap(double(ZA)) title({'Cameroon and Surrounding Region', aster.LayerTitle});

vmap0 = wmsfind('vmap0.tiles', 'SearchField', 'serverurl'); boundaries = refine(vmap0, 'country_02'); B = wmsread(boundaries, 'Latlim', latlim, ... 'Lonlim', lonlim, 'CellSize', cellSize, 'ImageFormat','image/png'); ZB = ZA; ZB(B(:,:,1) < 250) = min(ZA(:));

figure('Renderer','zbuffer') worldmap('Cameroon') demcmap(double(ZA)) geoshow(ZB, RA, 'DisplayType', 'texturemap') title({'Cameroon and Country Boundaries', aster.LayerTitle});

figure('Renderer','zbuffer') worldmap('Cameroon') imagesc((ZA)) geoshow(ZB, RA, 'DisplayType', 'texturemap') title({'Cameroon and Country Boundaries', aster.LayerTitle});

% figure % load korea % worlddatamap(map, refvec,'dem')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% S = shaperead('cntry02.shp','UseGeoCoords',true); % id = find(strcmpi('cameroon',{S.CNTRY_NAME})); % figure % latlim = S(19).BoundingBox(:,2); % lonlim = S(19).BoundingBox(:,1); % worlddatamap({'cntry02.shp'}, latlim, lonlim)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% import shapefile

%states = shaperead('CAM_boundaries', 'UseGeoCoords', true);

states = shaperead('CAM_outline', 'UseGeoCoords', true);

%states = shaperead('usastatelo', 'UseGeoCoords', true); lat = [states.Lat]; lon = [states.Lon]; [Z, R] = vec2mtx(lat, lon, 5, 'filled'); figure; worldmap(Z, R); geoshow(Z, R, 'DisplayType', 'texturemap') colormap(flag(3))

%%% TRIM ..... % maptrim(lat,lon)

% cameroon = shaperead('cntry98.shp',... % 'UseGeoCoords', true,... % 'Selector', {@(name)strcmpi('cameroon',name), 'Name'}); % inLat = cameroon.Lat; % inLon = cameroon.Lon;

% figure, % % imagesc(ZA) % hold on % plot(lon,lat,'.r')

[Jan]

@Matteo: Do you see, that your code is unreadable? If you see this, why don't you try to fix this?

Answer 1

[Cedric]

What you want to do is a zonal statistics, but I don't know how to do it properly in MATLAB (I have some sub-optimal solution though, see below). In ArcGIS, you would enable the extension Spatial Analyst and use a tool called "Zonal Statistics as Table". You would set the zonal feature class to be a polygon feature class of countries (or even select/extract the polygon(s) that correspond to relevant country/ies), set the zonal field to be a field of country codes, and use the raster of elevations as the raster to summarize. The output would provide you with statistics (sum, mean, std, min, max, count, range, etc) of the raster of elevation, per country code (structured as a table).

You can somehow achieve this in MATLAB by converting the polygon feature class of countries into a raster of country codes, with the same resolution as the raster of elevation. You also build a raster of pixels uniformly equal to 1 with the same resolution. Then you use ACCUMARRAY to obtain sums of elevations and sums of 1's for each country code. These allow you to compute an average elevation for each country. You can shorten the computation time by selecting only pixels from the raster of elevation that correspond to pixels with the relevant country code in the raster of country codes.

Assuming you had the raster of country codes available as a numeric array, the raster of elevations available as a numeric array.. if the code corresponding to the country you are interested in were 147, you could do something like

 mask = codesRaster == 147 ;
 elevations = elevationsRaster(mask) ;

You would then have all elevation values that correspond the country 147 and you could perform any count/stat/other based on them.