function [vx,vy,vv,ex,ey] = readicevelocity_0670(x,y,filepath)

% This function reads average ice sheet velocity from Greenland
% https://nsidc.org/data/nsidc-0670/versions/1
% 
% kristin poinar
% for glaciology class midterm project

%%%%%%%%%%%%%%%%%%%% INPUTS INFO %%%%%%%%%%%%%%%%%%%%
% x         desired locations in x, meters (north polar stereographic / EPSG:3413)
% y         desired locations in y, meters (north polar stereographic / EPSG:3413)
%
%%%%%%%%%%%%%%%%%%%% FILES NEEDED %%%%%%%%%%%%%%%%%%%%
% greenland_vel_mosaic250_vx_v1.tif
% greenland_vel_mosaic250_vy_v1.tif
% from the url in the header
% this is about 500 MB of data

%%%%%%%%%%%%%%%%%%%% OUTPUTS INFO %%%%%%%%%%%%%%%%%%%%
% vx        velocity in x direction, m/yr
% vy        velocity in y direction, m/yr
% vv        speed, m/yr

% Read from the TIF files
% filepath = './nsidc-0670/';
filestem = 'greenland_vel_mosaic250_';
        dx = 250;
        dy = 250;
        vels.x = -645125 : dx : (859875-dx);
        vels.y = -3370125: dy :(-640125-dy);
        vels.vx = flipud(imread(strcat(filepath,filestem,'vx_v1.tif')));
        vels.vy = flipud(imread(strcat(filepath,filestem,'vy_v1.tif')));
        vels.ex = flipud(imread(strcat(filepath,filestem,'ex_v1.tif')));
        vels.ey = flipud(imread(strcat(filepath,filestem,'ey_v1.tif')));

% Clip to x, y window we want
%         j = find(vels.x > min(x) & vels.x < max(x));
%         i = find(vels.y > min(y) & vels.y < max(y));
%         vels.vx = vels.vx(i,j);
%         vels.vy = vels.vy(i,j);
%         vels.x = vels.x(j);
%         vels.y = vels.y(i);
        
% Remove flagged pixels of no data (-2e9)
        temp = find(vels.vx < -1e9);
        vels.vx(temp) = NaN; 
        vels.vy(temp) = NaN; 
        clear dx dy temp

% Interpolate to the input grid
        vx = double(interp2(vels.x,vels.y,vels.vx,x,y));
        vy = double(interp2(vels.x,vels.y,vels.vy,x,y));
        ex = double(interp2(vels.x,vels.y,vels.ex,x,y));
        ey = double(interp2(vels.x,vels.y,vels.ey,x,y));
        
% Total velocity = speed
        vv = sqrt(vx.^2 + vy.^2);