%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% For use in optimization problem for genetic algorithms where you
% search for the highest elevation on a section of earth.
%
% By: Kevin Passino
% Date: 3/12/01
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

% For a specific country or region, go do the data base with
% in the US Dept. of Commerce, National Inst. of 
% Geophysical Data and get into the MGG-02 5min data
% To make it interesting let's look at a section of the earth
% in the Andean mountain range, particularly around Colombia
% South America.

% Download the data as space separated ascii data.  Save it
% on the computer and name the file country.dat
% then in matlab say "load country.dat -ascii"
% Now the data needs to be re-defined (the ordering
% does not conform to how matlab plots it.  So let
% country=flipud(country) to get it strait then do a
% "save country" and it will create a matlab data file
% called country.mat. To load that file

load colom.mat

% Good color plot, with edge of land, in color

% Define the patch of land on the earth for which we have
% topographical data (here, Colombia)

northlat=13; % Number of degrees north of the equator
southlat=-5; % Relative to north
westlong=-84; % Relative to Meridian of Greenwich
eastlong=-64; % Relative to Meridian of Greenwich

long=westlong:5/60:eastlong; % Define the vector of longitudinal
							% coordinates for plotting (loaded
							% 5 min data so there is a data
							% point each 5 min of a degree)
lat=southlat:5/60:northlat; % Define the vector of lateral
							% points at which we have data

% Check the definitions of the axes.

[lat_size,long_size]=size(colom); % Find the number of points 
if lat_size ~= size(lat), display('data error'), end
if long_size ~= size(long), display('data error'), end

% Next, plot the topographical data

figure(1)
clf
contour(long,lat,colom,20); colormap(jet)
xlabel('Degrees Longitudinal (- = west of Meridian of Greenwich)')
ylabel('Degrees North of Equator (- = south)')
title('Topographical map of Colombia (meters elevation)')
hold on
contour(long,lat,colom,[0 0],'k-')  % Adds the shorelines to the ocean
                             % and the Carribean by plotting the points
							 % where the contour is at 0 elevation
plot(long,0*ones(size(long)),'k--')  % Plots the equator

% Test the interpolation on the map so that elevations can 
% be computed anywhere on the map

interplong=westlong:3/60:eastlong; % Create 3 minute data
interplat=southlat:3/60:northlat;

interp_elev=interp2(long,lat,colom,interplong,interplat','linear');
figure(2)
clf
contour(interplong,interplat,interp_elev,20); colormap(jet)
xlabel('Degrees Longitudinal (- = west of Meridian of Greenwich)')
ylabel('Degrees North of Equator (- = south)')
title('Topographical map of Colombia (meters elevation)')
hold on
contour(long,lat,colom,[0 0],'k-')  % Adds the shorelines to the ocean
                             % and the Carribean by plotting the points
							 % where the contour is at 0 elevation
plot(long,0*ones(size(long)),'k--')  % Plots the equator


% Then to plot additional contour lines for specific elevations

%contour(colom,[1000 1000],'k--')
%contour(colom,[1800 1800],'k--')

%II=find(colom>=1000 && colom<=1800);
%III=find(colom>1800);

% For black and white printers

figure(2)
clf
contour(-colom,20); colormap(gray)
hold on
contour(-colom,[0 0],'k-')
xlabel('Degrees Longitudinal (- = west of Meridian of Greenwich)')
ylabel('Degrees North of Equator (- = south)')
title('Topographical map of Colombia (meters elevation)')


% For a 3d plot

figure(3)
clf
surf(colom)
shading interp
colormap(copper)
xlabel('Degrees Longitudinal (- = west of Meridian of Greenwich)')
ylabel('Degrees North of Equator (- = south)')
title('Topographical map of Colombia (meters elevation)')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%