Version8

clear all %clears and closes previously open figures
close all
clc

%Reading the .nc file
FileNC = 'o3_surface_20180701000000.nc' ;

%Printing information of the file in command window
ncinfo(FileNC)
%Longitude
long = ncread(FileNC,'lon');
%Length of longitude
nx = length(long) ;
%Latitude
lat = ncread(FileNC,'lat'); 
%Length of Latitude
ny = length(lat) ;
%Total number of time hours in a day
time = ncread(FileNC,'hour');

%Longitude and latitude for plotting. Converting it to double values
lon=double(long);
lat=double(lat);

%Selecting the data type
disp('Please select which climate model you would like to view: ');
type=input('1- chimere_ozone \n2- emep_ozone \n3- ensemble_ozone \n4- eurad_ozone \n5- lotoseuros_ozone \n6- match_ozone \n7- mocage_ozone \n8- silam_ozone \n');

%Data having all the data in it
if type==1
Data = ncread(FileNC,'chimere_ozone') ;
elseif type==2
Data = ncread(FileNC,'emep_ozone') ;    
elseif type==3
Data = ncread(FileNC,'ensemble_ozone') ;    
elseif type==4
Data = ncread(FileNC,'eurad_ozone') ;    
elseif type==5
Data = ncread(FileNC,'lotoseuros_ozone') ;    
elseif type==6
Data = ncread(FileNC,'match_ozone') ;    
elseif type==7
Data = ncread(FileNC,'mocage_ozone') ;    
elseif type==8
Data = ncread(FileNC,'silam_ozone') ;   
    
end

% user interface - allows the user to select an option to their needs
disp('Please select your desired format: ')
userInput = ('a - Standard Display (Trichromacy)\nb - Monochromacy Display \nc - Protanopia Display  \n');
option = input(userInput,'s');
%Creating meshgrid
[X,Y] = meshgrid(lat,lon);
figure (1)
%Trichromacy: Ozone Overlay 
if (option == 'a' || option == 'A') % this figure will display if a is selected
    
    for i = 1:24 %Loop for iterating
        clf; % figure is cleared to display following hour
        
    % Creating the map of EUROPE
    worldmap('Europe'); % Part of Earth to be displayed

    load coastlines % loading coastlines on map
    plotm(coastlat,coastlon); %plotting coast data 
    land = shaperead('landareas', 'UseGeoCoords', true);%Plotting Land
    geoshow(gca, land, 'FaceColor', [0.5 0.7 0.5]);
    lakes = shaperead('worldlakes', 'UseGeoCoords', true); %Plotting Lakes
    geoshow(lakes, 'FaceColor', 'blue'); 
    rivers = shaperead('worldrivers', 'UseGeoCoords', true); %Plotting rivers
    geoshow(rivers, 'Color', 'blue'); 
    cities = shaperead('worldcities', 'UseGeoCoords', true); % Plotting cities
    geoshow(cities, 'Marker', '.', 'Color', 'red');
   
    % Plots the data and sets transparency and edges
    surfm(X, Y, Data(:,:,i)', 'EdgeColor', 'none',...
        'FaceAlpha', 0.7); % edge colour outlines the edges, 'FaceAlpha', sets the transparency
    colormap default; % Normal vision mode
    colorbar % Displaying colour bar 
    title(sprintf("Standard View: Hour " + i + "\n" )); % title (displays hour and the type of view)
    pause(0.1); % time for each hour - pauses for 1 second
    end 

% Monochromacy: Ozone Overlay
elseif (option == 'b' || option == 'B') % this figure will display if c is selected
    for i = 1:24 %Loop for iterating
        clf; % figure is cleared to display following hour

   % Creating the map of EUROPE
    worldmap('Europe'); % Part of Earth to be displayed

    load coastlines % loading coastlines on map
    plotm(coastlat,coastlon); %plotting coast data 
    land = shaperead('landareas', 'UseGeoCoords', true);%Plotting Land
    geoshow(gca, land, 'FaceColor', [0.5 0.7 0.5]);
    lakes = shaperead('worldlakes', 'UseGeoCoords', true); %Plotting Lakes
    geoshow(lakes, 'FaceColor', 'blue'); 
    rivers = shaperead('worldrivers', 'UseGeoCoords', true); %Plotting rivers
    geoshow(rivers, 'Color', 'blue'); 
    cities = shaperead('worldcities', 'UseGeoCoords', true); % Plotting cities
    geoshow(cities, 'Marker', '.', 'Color', 'red');

    % Plotting the data
    surfm(X, Y, Data(:,:,i), 'EdgeColor', 'none',...
        'FaceAlpha', 0.7) % edge colour outlines the edges, 'FaceAlpha', sets the transparency
    colormap gray; %ColorBlindness mode for monachromatic vision
    colorbar % Displaying colour bar 
    title(sprintf("Monochromacy View: Hour " + i + "\n" )) % title
    pause(0.1) % time for each hour
    end 
% Protanopia: Ozone Overlay
elseif (option == 'c' || option == 'C') % this figure will display if c is selected
    for i = 1:24 %Loop for iterating 
        clf; % figure is cleared to display following hour

    % Creating the map of EUROPE
    worldmap('Europe'); % Part of Earth to be displayed

    load coastlines % loading coastlines on map
    plotm(coastlat,coastlon); %plotting coast data 
    land = shaperead('landareas', 'UseGeoCoords', true);%Plotting Land
    geoshow(gca, land, 'FaceColor', [0.5 0.7 0.5]);
    lakes = shaperead('worldlakes', 'UseGeoCoords', true); %Plotting Lakes
    geoshow(lakes, 'FaceColor', 'blue'); 
    rivers = shaperead('worldrivers', 'UseGeoCoords', true); %Plotting rivers
    geoshow(rivers, 'Color', 'blue'); 
    cities = shaperead('worldcities', 'UseGeoCoords', true); % Plotting cities
    geoshow(cities, 'Marker', '.', 'Color', 'red');

    % Plot the data
    surfm(X, Y, Data(:,:,i), 'EdgeColor', 'none',...
        'FaceAlpha', 0.7) % edge colour outlines the edges, 'FaceAlpha', sets the transparency
    colormap Jet ; % ColorBlindness mode for protanopes
    colorbar % diaplaying colorbar
    title(sprintf("Protanopia View: Hour " + i + "\n" )) % title
    pause(0.1) % time for each hour
    end 
   
else 
    %Error Message if no option selected
    disp('Sorry that is not an option from above! Please run the program again')
    
end