MATLAB: How to filter the required coordinates which consist with same elevation

filterlat/lon coordinates

Hi there,

I have a netCDF dataset which consists of longitudes (1xN), latitudes (Mx1) and different elevations (MxN). I need to filter the lat/lon coordinates which are having zero (0 m) elevation and save those coordinates in a seperate text file.

I have used following line of codes to extract the netCDF data, but struggling the filtering part.

filename = 'elevations.nc';
elev = ncread(filename, 'elevation')';
lon  = ncread(filename, 'lon')';
lat  = ncread(filename, 'lat');

I would greatly appreciate anyone can help me to solve this problem.

Thank you.

Best Answer

filename = 'elevations.nc';
elev = ncread(filename, 'elevation')';
lon  = ncread(filename, 'lon')';
lat  = ncread(filename, 'lat');
[c,h] = contour(lon,lat,elev,'showtext','on','LevelList',0) ;

YOu have your o contour coordinates in c.

Related Solutions

MATLAB: Visualizing Geoid Height with given longitude, latitude and geoid heights

First, create a 2D map of the heights were the rows are latitude and the columns are longitude.

You could use the matlab unique() on the list of latitudes and on the list of longitudes to create 1D arrays to define the row and column values. Then use the meshgrid command to make a 2D map of latitude and a 2D map of longitude. If the data are nicely ordered you can just use the reshape command to make the 1D array of heights into a 2D array. If it is not, then loop through lat and lon of the 2D array (each pixel) and search the Nx2 array accordingly doing a nearest neighbor interpolation into the 2D array of heights.

Most likely the data are nicely ordered so that reshape will work to create the 2D array. You might need to transpose the array. If each record first loops through a series of latitudes for one longitude, then increments longitude and loops through the same set of latitudes, then the height data. Below is some script code that creates fake topo data to illustrate:

nlong = 50;

nlat = 25;

[ICOL,IROW] = meshgrid(1:nlong,1:nlat);

longitude = zeros(nlat,nlong);

latitude = zeros(nlat,nlong);

height = zeros(nlat,nlong);

for ilong = 1:nlong

for ilat = 1:nlat

longitude(ilat,ilong) = 0.5*ilong - 70;

latitude(ilat,ilong) = 0.5*ilat + 40;

height(ilat,ilong) = 2i*pi*rand(1,1);

end

% making better fake data to illustrate

height = abs(ifft2(ifftshift((hypot(IROW,ICOL)+1e-5).^(-2.1).*exp(height))));

% shows ordering of longitude outer cycle and latitude inner cycle

% this is very likely how your data are organized

X = reshape(longitude,nlat*nlong,1);

Y = reshape(latitude,nlat*nlong,1);

Z = reshape(height,nlat*nlong,1);

for ipt = 1:nlat*nlong

fprintf(1,'%f %f %f\n',X(ipt),Y(ipt),Z(ipt));

end

% if this is how your data are arranged, then to create the 2D map:

latitude1D = unique(Y);

longitude1D = unique(X);

height2D = reshape(Z,nlat,nlong);

latitude2D = reshape(Z,nlat,nlong);

longitude2D = reshape(Z,nlat,nlong);

% if the reverse is true, that latitude is the outer cycle, then you need

% to transpose the result

% height2D = height2D';

%% referene matrix

% Reference matrix for mapping geoid heights to lat/lon on globe.

R1 = makerefmat('RasterSize',size(height2D), ...

'Latlim', [min(latitude(:)) max(latitude(:))], ...

'Lonlim', [min(longitude(:)) max(longitude(:))] );

%% visualization 1: DTED map with contour lines

figure();

hold on;

mapshow(zeros(size(height2D)),R1,'CData',height2D,'DisplayType','surface');

cW = mapshow(height2D,R1,'DisplayType','contour','linecolor','black', ...

'showtext','on');

demcmap(height2D);

axis equal

axis tight

%% visualization 2: with coastline

coast = load('coast');

% Create 2-D plot using |meshm|

h2D = figure;

set(h2D,'Position',[20 75 700 600],'Toolbar','figure');

gridDegInc = 0.5;

ast2DGeoidPlot(R1,height2D,coast,gridDegInc)

%% another visualization using the aerospace toolbox

figure

gcm = usamap([min(latitude1D) max(latitude1D)], [min(longitude1D) max(longitude1D)]);

hold on;

geoshow(gcm,height2D,R1,'DisplayType','texturemap');

demcmap(height2D);

daspectm('m',1);

geoshow(latitude2D,longitude2D,height2D,'DisplayType','contour', ...

'LineColor','black');

%% another visualization

figure

gcm = usamap([min(latitude1D) max(latitude1D)], [min(longitude1D) max(longitude1D)]);

hold on;

geoshow(gcm,height2D,R1,'Zdata',height2D, 'CData',height2D,'DisplayType','surface')

demcmap(height2D)

daspectm('m',1)

view(3)

MATLAB: How to extract the data from .nc file based on longitude and latitude,time

The nc* functions and the netcdf.* functions represent two independent ways to access a NetCDF file. I would recommend sticking to nc* functions until you get comfortable enough to directly access the NetCDF library interface using the netcdf.* functions (it uses C style indexing and dimension ordering which is different from what nc* functions report and use).

Instead of

 precJanuary = netcdf.getVar(ncid,2,[0 0 0],[128 64 1]);

try

 precJanuary = ncread(filename, varname,[1 1 1],[128 64 1]);

Best Answer

Related Solutions

MATLAB: Visualizing Geoid Height with given longitude, latitude and geoid heights

MATLAB: How to extract the data from .nc file based on longitude and latitude,time

Related Question