%Getting the velocity in PIG vicinity for the ExpDraw

% Load Velocities
% http://nsidc.org/data/nsidc-0484.html
nsidc_vel='../Data/Antarctica_ice_velocity.nc'; 	

% Get necessary data to build up the velocity grid
xmin = ncreadatt(nsidc_vel,'/','xmin');
xmin = strtrim(xmin);  % this is a string, and we need to recover the double value
xmin = xmin(1:end-2);  % get rid of the unit
xmin = str2num(xmin);  % convert to double

ymax = ncreadatt(nsidc_vel,'/','ymax'); 
ymax = strtrim(ymax);  
ymax = ymax(1:end-2);  
ymax = str2num(ymax); 
	
nx = ncreadatt(nsidc_vel,'/','nx');
ny = ncreadatt(nsidc_vel,'/','ny');

spacing = ncreadatt(nsidc_vel,'/','spacing'); 
spacing = strtrim(spacing);
spacing = spacing(1:end-2);  
spacing = str2num(spacing); 


x=xmin+(0:1:nx-1)'*spacing; 
x=double(x);
y=(ymax)-(0:1:ny-1)'*spacing; 
y=double(y);

posx=find(x<=-12.0e5);
id1x=find(x>=-18.0e5,1);
id2x=posx(end);

posy=find(y>=-4.0e5);
id1y=find(y<=1.0e5,1);
id2y=posy(end);

%Get velocity subset
vx = double(ncread(nsidc_vel,'vx'));
vx_obs = vx(id1x:id2x,id1y:id2y);
vx_obs = flipud(vx_obs');

vy = double(ncread(nsidc_vel,'vy'));
vy_obs = vy (id1x:id2x,id1y:id2y);
vy_obs = flipud(vy_obs');

xred=x(id1x:id2x);
yred=y(id1y:id2y);
vel_obs=sqrt(vx_obs.^2.+vy_obs.^2.);
imagesc(xred,yred,vel_obs)