% load your data
%x = vector of longitude
%y = vector of latitude
%u = u-velocity (dimensions: longitude x latitude)
%v = v-velocity (dimensions: longitude x latitude)

pkg load netcdf

t = cputime;

filename = 'sscurrentvel_annual.nc'; 
delete(filename);

%x=[-170:1:170]';
%y=[-80:1:80]';

%u=rand(size(x,1),size(y,1));
%v=rand(size(x,1),size(y,1));

fprintf('Reading values from files\n');
drawnow();

#TODO: read file and assign NA to _
%zvelocity=dlmread("zvelocity.csv");
%mvelocity=dlmread("mvelocity.csv");

fileID = fopen("zvelocity.csv");
zvelocity2 = textscan(fileID,"%f %f %f %s",'delimiter',' ','EmptyValue',NaN);
fclose(fileID);
zvelocity2{1,4}=strrep(zvelocity2{1,4},"_","0");
latRaw=(zvelocity2{1,2}(:));
longRaw=(zvelocity2{1,3}(:));
qRaw=str2double(zvelocity2{1,4}(:,1)');

uall=[];
uall(:,1)=latRaw;
uall(:,2)=longRaw;
uall(:,3)=qRaw;

fileID = fopen("mvelocity.csv");
mvelocity2 = textscan(fileID,"%f %f %f %s",'delimiter',' ','EmptyValue',NaN);
fclose(fileID);
mvelocity2{1,4}=strrep(mvelocity2{1,4},"_","0");
latRaw=(mvelocity2{1,2}(:));
longRaw=(mvelocity2{1,3}(:));
qRaw=str2double(mvelocity2{1,4}(:,1)');

vall=[];
vall(:,1)=latRaw;
vall(:,2)=longRaw;
vall(:,3)=qRaw;


%uall=zvelocity(:,2:4);
%vall=mvelocity(:,2:4);

y=sort(unique(uall(:,1)));
x=sort(unique(uall(:,2))-180);

lat=sort(unique(uall(:,1)));
long=sort(unique(uall(:,2)));

u = zeros(length(y),length(x));
v = zeros(length(y),length(x));
u_non_proj = zeros(length(y),length(x));
v_non_proj = zeros(length(y),length(x));


fprintf('Assigning values to variables\n');
drawnow();

np=length(long)*length(lat);
nl = length(long);
nlh = round(nl/2);

for i = 1:length(long)
    latuvaluesI=find(uall(:,2)==long(i));
    latvvaluesI=find(vall(:,2)==long(i));

    latuvalues=uall(latuvaluesI,1:end);
    latvvalues=vall(latvvaluesI,1:end);
    
    latuvalues=sortrows(latuvalues,1);
    latvvalues=sortrows(latvvalues,1);
    uvalues = latuvalues(:,3);
    %uvalues(uvalues==0) = NA;
    vvalues = latvvalues(:,3);
    %vvalues(vvalues==0) = NA;


    u_non_proj(:,i)=uvalues;
    v_non_proj(:,i)=vvalues;
 
    %u(:,i)=latuvalues(:,3);
    %v(:,i)=latvvalues(:,3);
end

for i = 1:length(lat)
for j = 1:nl
	k = mod((j+(nlh)),nl);
	if (k>nl)
		k=nl;
	end
	if (k<1)
		k=1;
	end
	u(i,j)=u_non_proj(i,k);
	v(i,j)=v_non_proj(i,k);
end
end

fprintf('Finished assigning values to variables\n');
drawnow();
disp('Elapsed: ');
e = cputime-t

nccreate(filename,'lon','Dimensions',{'lon',size(u,2)},'Format','classic');
ncwriteatt(filename,'lon','standard_name','longitude');
ncwriteatt(filename,'lon','units','degree_east');

nccreate(filename,'lat','Dimensions',{'lat',size(u,1)});
ncwriteatt(filename,'lat','standard_name','latitude');
ncwriteatt(filename,'lat','units','degree_north');

nccreate(filename,'u','Dimensions',{'lon','lat'});
ncwriteatt(filename,'u','standard_name','northward_sea_water_velocity');
ncwriteatt(filename,'u','units','m s-1');

nccreate(filename,'v','Dimensions',{'lon','lat'});
ncwriteatt(filename,'v','standard_name','eastward_sea_water_velocity');
ncwriteatt(filename,'v','units','m s-1');

ncwriteatt(filename,'/','Conventions','CF-1.5');

#note that the matrices should be trasposed
ncwrite(filename,'lon',x);
ncwrite(filename,'lat',y);
ncwrite(filename,'u',u'); 
ncwrite(filename,'v',v');