From: Paul Alexander on
I would appreciate some help in what seems to be a simple problem that I'm having a lot of trouble with.
I have longitude and latitude in gridded files. ( Along with a bevy of other gridded files like land/sea masks, bathymetry, etc.) I need to be able to read in a long x lat file that defines a boundary or region of the ocean to be simply identified by a number. The end result of this should be a grid file of the same size as the lat and long grids but instead of locations I need it to just list numbers. The real problem is that I need to be able to compare where these boundary lines are in the real world to where the computer represents them with the limited data of a grid.
For example if I put in a boundary file that encompasses some polygon of space in the real world I need to be able to represent that in Matlab with a matrix that corresponds to pre-existing lat and long grid files. Any ideas?
If I haven't been clear enough please advise and I'll rewrite this to offer a better explanation of my problem.
From: Paul Alexander on
Here is the code I've written. The index variable "same" doesn't return any values. It's kind of convoluted but what I'm trying to do from lines 74 to 78 is find where my boundary line is either A) less than both the long and lat values at increment 'j' and greater than than both the lat and long values at increment 'j+1' or B) greater than both the long and lat values at increment 'j' and less than both the lat and long values at increment 'j+1'. I had hoped that would yield the places where the boundary line is in between the real world locations but it returns nothing. Any thoughts?

%-------------------------------------------------------------------------%
% Get input to determine how many regions are needed. Boundary line data %
% will need to be in seperate files. %
%-------------------------------------------------------------------------%
iter=input('How many regions? ');

%-------------------------------------------------------------------------%
% Get the boundary line files and evaluate. %
%-------------------------------------------------------------------------%
for k=1:iter;
bdyname= input(['Enter boundary line (lat,lon) file #'...
num2str(k) ' '],'s');
[templat,templon]=rcoastline(bdyname);
%-------------------------------------------------------------------------%
% Convert longitudes >180 to negative longitudes. %
%-------------------------------------------------------------------------%
for j=1:size(templon);
if templon(j,:)>=180 && ~isnan(templon(j,:));
templon(j,:)=templon(j,:)-360;
end
end
clear j
%-------------------------------------------------------------------------%
% Assign longitude and latitude values to variables. %
%-------------------------------------------------------------------------%
eval(['bdylon' num2str(k) ' =templon;']);
eval(['bdylat' num2str(k) ' =templat;']);
clear bdyname templat templon
end
clear k
%-------------------------------------------------------------------------%
% Use 'read_mask.m' to get variables: spherical, rlon, rlat, bath & mask. %
%-------------------------------------------------------------------------%
grdname=input('Please enter the GRID netCDF file name: ','s');
[spherical,rlon,rlat,bath,mask]=read_mask(grdname);
%-------------------------------------------------------------------------%
% Create a grid the same size as 'mask' but for the use of region numbers.%
%-------------------------------------------------------------------------%
regions=zeros(size(mask));
%-------------------------------------------------------------------------%
% Assign temporary variables (rlonv and rlatv) a vectorized format of rlon%
% and rlat for comparison against the boundary line values. %
%-------------------------------------------------------------------------%
rlonv=reshape(rlon, numel(rlon),1);
rlatv=reshape(rlat, numel(rlat),1);
%-------------------------------------------------------------------------%
% In order to compare boundary lines and rlon and rlat we need to truncate%
% after the 2nd decimal (because otherwise the points would be too precise%
% to ever actually be the same). So I matrix multiplied by 100 and casted %
% the resulting matrix as an int16. %
%-------------------------------------------------------------------------%
rlatv=rlatv*100;
rlatv=int16(rlatv);
rlonv=rlonv*100;
rlonv=int16(rlonv);
%-------------------------------------------------------------------------%
% Create a temporary variable that allows you to assign value and index %
% but will still change with every iteration. Also truncate the variable. %
%-------------------------------------------------------------------------%
for k=1:iter;
templon=eval(['bdylon' num2str(k)]);
templat=eval(['bdylat' num2str(k)]);
templon=templon*100;
templon=int16(templon);
templat=templat*100;
templat=int16(templat);
%-------------------------------------------------------------------------%
% Compare every value in the temporary longitude and latitude to every %
% value in the vectorized rho-longitude and rho-latitude. Index where %
% the templon and templat are inbetween the rho-long and rho-lat values. %
%-------------------------------------------------------------------------%
for i=1:size(templon)
for j=1:(1-size(rlonv));
same=(((templon(i,:)<=rlonv(j,:)&&templon(i,:)>=rlonv(j+1,:))...
||(templon(i,:)>=rlonv(j,:)&&templon(i,:)<=rlonv(j+1,:)))...
&&((templat(i,:)<=rlatv(j,:)&&templat(i,:)>=rlatv(j+1,:))...
||(templat(i,:)>=rlatv(j,:)&&templat(i,:)<=rlatv(j+1,:))));
%-------------------------------------------------------------------------%
% Assign region numbers to indexed locations in a new mask grid. %
%-------------------------------------------------------------------------%
regions(same)=k;
end
clear j same
end
clear i
end
clear k j i templon templat same