function OK = convertCoordinates_test(varargin)
%CONVERTCOORDINATES2_TEST   test convertCoordinates
% 
% Several tests to check the coordinate conversion routines. 
% The EPSG database s loaded only once for speed. Not all tests have a hard
% criterium, the first tests passes if only 75% percent of operations
% succeed. 
% 
%See also: CONVERTCOORDINATES, CONVERTCOORDINATES_TEST

MTestCategory.DataAccess;

% Test is shortened to run in 1 minute
% if TeamCity.running
%     TeamCity.ignore('Test takes very long');
%     OK = 1;
%     return;
% end

   if nargin==1
   EPSG = varargin{1};
   else
   EPSG = load('EPSG');
   end

   EPSG        = load('EPSG');
   OK(1)       = test1(EPSG);
   OK(2)       = test2(EPSG);
   
   OK(3)       = convertCoordinates_test_netherlands(EPSG);
   OK(4)       = convertCoordinates_test_germany    (EPSG);
  %OK(5)       = convertCoordinates_test_uk         (EPSG);
   OK          = all(OK);

function OK = test1(EPSG)
% This test converts from all systems to GEO-WGS84 and back and it 
% checks than if the output compares well with the input coordinates.
% Results are printed to output file convertCoordinate_errorMessages.txt
% The succes criteria is that this goes well on at least 75% of the
% conversions

% load coordinate systems
codes       = EPSG.coordinate_reference_system.coord_ref_sys_code;
kinds       = EPSG.coordinate_reference_system.coord_ref_sys_kind;
deprecated  = EPSG.coordinate_reference_system.deprecated;
names       = EPSG.coordinate_reference_system.coord_ref_sys_name;

% only check non-deprecated geographic 2D and projected sytems, ignore others
codes = codes(ismember(kinds,{'geographic 2D', 'projected'}) & strcmpi(deprecated,'FALSE'));
names = names(ismember(kinds,{'geographic 2D', 'projected'}) & strcmpi(deprecated,'FALSE'));

geo_code = 4326; % 4326 = WGS 84

% create file for printing error messages
pat = fileparts(which('convertCoordinates_test'));
fid = fopen([pat filesep 'convertCoordinate_errorMessages' version '.txt'],'w');

% preacclocate testresult
testresult = nan(size(codes));

% start loop over all systems
for s = 1:length(codes)
    fprintf(fid,'Code :%-9.0f| Name:%-68s| ',codes(s), names{s});
    try
        % try conversion at the center and the limits of the
        % area_of_use (area where coordinate system should be valid)
        
        ind           = find(EPSG.coordinate_reference_system.coord_ref_sys_code == codes(s),1,'first');
        ind           = find(EPSG.area.area_code == EPSG.coordinate_reference_system.area_of_use_code(ind),1,'first');
        
        [lat,lon]     = meshgrid(...
            linspace(EPSG.area.area_north_bound_lat(ind),EPSG.area.area_south_bound_lat(ind),3),...
            linspace(EPSG.area.area_east_bound_lon(ind) ,EPSG.area.area_west_bound_lon(ind) ,3));
        lat = lat([5 1 3 7 9]); % center and corners
        lon = lon([5 1 3 7 9]); % center and corners
        
        [x   ,y   , OPT]   = convertCoordinates(lon ,lat ,EPSG,'CS1.code', geo_code,'CS2.code', codes(s),...
            'CS1.UoM.name','degree'); % with 'CS1.UoM.name','degree' the input is forced to be in degrees
        
        [lon2,lat2]        = convertCoordinates(x   ,y   ,EPSG,'CS2.code', geo_code,'CS1.code',codes(s),...
            'CS2.UoM.name','degree',...       % make sure the output is in degree, and the input
            'CS1.UoM.name',OPT.CS2.UoM.name); % is in the same system as the output of the previous conversion
        
        testresult(s) = all(abs(lon - lon2) < 0.0001 & abs(lat -lat2) < 0.0001);
        if testresult(s)
            fprintf(fid,'Test result is OK!\n');
        else
            % check if correct at the center
            if abs(lon(1) - lon2(1)) < 0.0001 && lat(1) -lat2(1) < 0.0001
                fprintf(fid,'Test result is OK in center, but not in corners: diff(lon) =%s | diff(lat) =%s \n',sprintf('% 10.4f',lon(:)-lon2(:)), sprintf('% 10.4f',lat(:)-lat2(:)));
            else
                fprintf(fid,'Test result is not OK: diff(lon) =%s | diff(lat) =%s \n',sprintf('% 10.4f',lon(:)-lon2(:)), sprintf('% 10.4f',lat(:)-lat2(:)));
            end
        end
    catch
        testresult(s) = nan;
        fprintf(fid,'%s\n',strrep(lasterr,char(10),' '));
    end
end

fprintf(fid,'\nTotal score: %-4.1f%% passed this simple test',sum(testresult==1)/numel(testresult)*100);

fclose(fid);

OK = sum(testresult==1)/numel(testresult)*100 > 75;

function OK = test2(EPSG)
% test to see if unit conversions work

geo_code = 4326; % 4326 = WGS 84

lon = [-116.5950 -115.0000 -118.1900 -115.0000 -118.1900];
lat = [  38.5000   41.0000   41.0000   36.0000   36.0000];

[x1,y1]   = convertCoordinates(lon ,lat ,EPSG,'CS1.code', geo_code,'CS2.name', 'NAD83(NSRS2007) / Nevada Central');

[x2,y2]   = convertCoordinates(lon ,lat ,EPSG,'CS1.code', geo_code,'CS2.name', 'NAD83(NSRS2007) / Nevada Central (ft US)', 'CS2.UoM.name','metre'); % with 'CS2.UoM.name', the output is forced to be in metres

OK = all(x1-x2 < 0.001) && all(y1-y2 < 0.001);
