function varargout = obs_file_in_polygon(varargin);
%obs_file_in_polygon - construct an obervation file (*_obs.xyn) from all locations within a polygon
%
% Creates a D-Flow FM observation file, based on all flownodes within 
% (and on) a predefined polygon/region.
%
%    <output> = dflowfm.obs_file_in_polygon(<1>,<2>,<3>);
%
%    <1>      = D-Flow FM *_map.nc file (optional)
%    <2>      = Polygon *.pol or landboundary *.ldb file (optional)
%    <3>      = Observation filename *_obs.xyn (optional)
%    <output> = [Mx2] matrix with locations of observations (optional)
%
%   Note that ALL input variables are optional, when simply calling
%   dflowfm.obs_file_in_polygon the user is asked to supply the *_map.nc
%   file & polygon/landboundary file (see examples 1 and 2)
%
% Input variable specification:
%
%  <D-Flow FM *_map.nc file>
%    Can either be a string containing the location of the *_map.nc file
%    (both with or without the folder)or a structure generated by
%    dflowfm.readNet (note that only *_map.nc files are supported, as grid
%    files (*_net.nc) do not contain flownode locations)
%
%  <Polygon *.pol or landboundary *.ldb file>
%    Can either be a string containing the location of a *.pol or *.ldb
%    file (both with or without the folder), or an [Mx2] or [2xM] matrix
%    (generated by e.g. the functions landboundary, LDBTool or manually)
%
%  <Observation filename *_obs.xyn>
%    String containing the name of the observation file (output), must end 
%    with _obs.xyn, else, the default name is used instead (which is 
%    output_obs.xyn, in the current directory). A directory can also be
%    specified along with the name, e.g. 'D:\output\example_obs.xyn'
%
% Output variable specification:
%
%  <Output>
%    [Mx2] matrix, containing all locations of observation points
%
%
%  Example 1:
%    obs_locs = dflowfm.obs_file_in_polygon; *
%              * (User will be asked for all files)
%
%  Example 2:
%    dflowfm.obs_file_in_polygon('D:/dflowfm/output_map.nc'); *
%              * (User will be asked for a polygon/landboundary)
%
%  Example 3:
%    dflowfm.obs_file_in_polygon('D:/dflowfm/output_map.nc','D:/dflowfm/ldb_file.ldb','filename_obs.xyn');
%
%  Example 4:
%    grid_data = dflowfm.readNet('D:/dflowfm/output_map.nc');
%    poly_data = landboundary('read','D:/dflowfm/polygon_file.pol');
%    obs_locs  = dflowfm.obs_file_in_polygon(grid_data,poly_data);
%    dflowfm.plotNet(grid_data); hold on;
%    plot(obs_locs(:,1),obs_locs(:,2),'.');
%
%  Example 5:
%    grid_data = dflowfm.readNet('D:/dflowfm/output_map.nc');
%    dflowfm.obs_file_in_polygon(grid_data,[1 2; 1 3; 3 3; 2 1],'D:\output\filename_obs.xyn');
%
%
%  Additional note (advanced users):
%    If the user allows for observation points to be placed on netnode
%    locations (so interpolation of the modelling results towards these
%    points will be made), simply generate a grid structure from the
%    *net.nc grid file using: grd = dflowfm.readNet('grid_net.nc') and then
%    state: grd.cen = grd.cor to fool the function in advance. Note that
%    this option removes the need for a *_map.nc file, though the locations
%    of observation points are no longer on flownodes (might not be
%    preferable).
%
% See also: dflowfm, dflowfm.readNet, inpolygon, landboundary, LDBTool, delft3d_obs_file_in_polygon

%   --------------------------------------------------------------------
%   Copyright (C) 2013 Deltares
%       Freek Scheel
%
%       <freek.scheel@deltares.nl>;
%
%       Deltares
%       P.O. Box 177
%       2600 MH Delft
%       The Netherlands
%
%   This library is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this library.  If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

if nargout>1
    error('Too many output arguments, use either none or one');
end

if size(varargin,2)>3
    error('Too much input defined');
elseif isempty(varargin)
    uiwait(msgbox('No D-Flow FM map-file and polygon/landboundary defined, please select these now...','Map-file and polygon missing','modal'))
    [map_name, map_path] = uigetfile('*_map.nc', 'Select a D-Flow FM map output file (*_map.nc)');
    if size(map_name==0,2)==1
        if map_name==0
            error('Aborted by user, since no *_map.nc file was selected');
        end
    end
    [pol_name, pol_path] = uigetfile({'*.pol;*.ldb'}, 'Select a polygon (*.pol) or landboundary (*.ldb) file');
    if size(pol_name==0,2)==1
        if pol_name==0
            error('Aborted by user, since no *.pol or *.ldb file was selected');
        end
    end
    
    tic
    
    map_file_loc = [map_path map_name];
    map_file = dflowfm.readNet(map_file_loc);
    
    pol_file_loc = [pol_path pol_name];
    pol_file = landboundary('read',pol_file_loc);
    
    obs_file_loc = [pwd filesep 'output_obs.xyn'];
elseif size(varargin,2)==1
    uiwait(msgbox('No polygon/landboundary file defined, please select this now...','Polygon/landboundary missing','modal'))
    [pol_name, pol_path] = uigetfile({'*.pol;*.ldb'}, 'Select a polygon (*.pol) or landboundary (*.ldb) file');
    if size(pol_name==0,2)==1
        if pol_name==0
            error('Aborted by user, since no *.pol or *.ldb file was selected');
        end
    end
    
    tic
    
    if ischar(varargin{1})==1
        if isempty(strfind(varargin{1},filesep))==1;
            map_file_loc = [pwd filesep varargin{1}];
            map_file = dflowfm.readNet(map_file_loc);
        elseif ~isempty(strfind(varargin{1},filesep))==1;
            map_file_loc = [varargin{1}];
            map_file = dflowfm.readNet(map_file_loc);
        end
    elseif isstruct(varargin{1})==1
        map_file_loc = 0;
        if isfield(varargin{1},'cen')==1
            map_file = varargin{1};
        elseif isfield(varargin{1},'cor')==1
            error('It appears that you selected a D-Flow FM grid file\nPlease select a D-Flow FM map output file instead','error');
        else
            error('Not a valid file format choosen');
        end
    else
        help obs_file_in_polygon
        error('Please specify your D-Flow FM map-file function input format according to the options given by the help above');
    end
    
    pol_file_loc = [pol_path pol_name];
    pol_file = landboundary('read',pol_file_loc);
    
    obs_file_loc = [pwd filesep 'output_obs.xyn'];
elseif size(varargin,2)>1
    
    tic
    
    if ischar(varargin{1})==1
        if isempty(strfind(varargin{1},filesep))==1;
            map_file_loc = [pwd filesep varargin{1}];
            map_file = dflowfm.readNet(map_file_loc);
        else ~isempty(strfind(varargin{1},filesep))==1;
            map_file_loc = [varargin{1}];
            map_file = dflowfm.readNet(map_file_loc);
        end
    elseif isstruct(varargin{1})==1
        map_file_loc = 0;
        if isfield(varargin{1},'cen')==1
            map_file = varargin{1};
        elseif isfield(varargin{1},'cor')==1
            error({'It appears that you selected a D-Flow FM grid file';'Please select a D-Flow FM map output file instead'});
        else
            error('Not a valid file format choosen');
        end
    else
        help obs_file_in_polygon
        error('Please specify your D-Flow FM map-file function input format according to the options given by the help above');
    end
    
    if ischar(varargin{2})==1
        if isempty(strfind(varargin{2},filesep));
            pol_file_loc = [pwd filesep varargin{2}];
            pol_file = landboundary('read',pol_file_loc);
        else ~isempty(strfind(varargin{2},filesep));
            pol_file_loc = [varargin{2}];
            pol_file = landboundary('read',pol_file_loc);
        end
    elseif ismatrix(varargin{2})==1
        pol_file_loc = 0;
        if size(varargin{2},2)==2
            pol_file = varargin{2};
        elseif size(varargin{2},1)==2
            pol_file = varargin{2}';
        else
            error('Inconsistent polygon/landboundary format, only [X,Y] formats are supported');
        end
    else
        help obs_file_in_polygon
        error('Please specify your polygon/landboundary function input format according to the options given by the help above');
    end
    obs_file_loc = [pwd filesep 'output_obs.xyn'];
end

if size(varargin,2)==3
    if ischar(varargin{3})
        if size(varargin{3},2)>8
            if strcmp(varargin{3}(end-7:end),'_obs.xyn')==1
                if isempty(strfind(varargin{3},filesep));
                    obs_file_loc = [pwd filesep varargin{3}];
                elseif ~isempty(strfind(varargin{3},filesep));
                    obs_file_loc = [varargin{3}];
                end
            else
                warning('Note that the obs file should end with _obs.xyn, using default filename instead (output_obs.xyn)');
                obs_file_loc = [pwd filesep 'output_obs.xyn'];
            end
        else
            warning('Please specify an output name ending with _obs.xyn, default filename is used instead (output_obs.xyn)');
            obs_file_loc = [pwd filesep 'output_obs.xyn'];
        end
    else
        warning('Please specify the output file name as characters, default filename is used instead (output_obs.xyn)');
        obs_file_loc = [pwd filesep 'output_obs.xyn'];
    end
end

[in_p on_p] = inpolygon(map_file.cen.x,map_file.cen.y,pol_file(:,1),pol_file(:,2));
in_poly_final = logical(in_p + on_p);

if size(in_poly_final,2)>99999
    error('Too many points found, over 100.000!');
end

obs_X = map_file.cen.x(in_poly_final); obs_Y = map_file.cen.y(in_poly_final);

obs_file = fopen(obs_file_loc,'w');
if max(max([obs_X obs_Y]))<180
    for ii = 1:size(obs_X,2)
        fprintf(obs_file,['%12.8f %12.8f obs_no_%05.0f\n'],obs_X(ii),obs_Y(ii),ii);
    end
else
    for ii = 1:size(obs_X,2)
        fprintf(obs_file,['%12.3f %12.3f obs_no_%05.0f\n'],obs_X(ii),obs_Y(ii),ii);
    end
end
fclose(obs_file);

obs_file_loc_fileseps = strfind(obs_file_loc,filesep);

if nargout==1
    varargout{1} = [obs_X' obs_Y'];
end

if size(obs_X,2)>1
    disp(' ');
    disp(['Succesfully generated observation file ' obs_file_loc(obs_file_loc_fileseps(end)+1:end) ' with ' num2str(size(obs_X,2)) ' observation points in ' num2str(round(toc)) ' seconds..']);
    disp(['File can be found in folder ''' obs_file_loc(1:obs_file_loc_fileseps(end)) '''']);
else
    disp(' ');
    disp(['Succesfully generated observation file ' obs_file_loc(obs_file_loc_fileseps(end)+1:end) ' with ' num2str(size(obs_X,2)) ' observation point in ' num2str(round(toc)) ' seconds..']);
    disp(['File can be found in folder ''' obs_file_loc(1:obs_file_loc_fileseps(end)) '''']);    
end











