|
|
trimming a matrix (raster grid) with a shapefile ....
Posted:
Mar 10, 2013 11:45 AM
|
|
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
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')
|
|
