% TG_LIDAR_COMP:  Reads in radar and lidar DEMs, compares slopes along 
% flowline.

% Read in files.  Lidar DEM first:

fid2                  = fopen('./DEM_lidar_utm57s.dat','r');
 
breaktest              = 0;
i                      = 0;

new_line21             = fgets(fid2);
[dummy21,ncols2]       = strread(new_line21,'%s %f','whitespace',' ');
new_line22             = fgets(fid2);
[dummy22,nrows2]       = strread(new_line22,'%s %f','whitespace',' ');
new_line23             = fgets(fid2);
[dummy23,xllcorner2]   = strread(new_line23,'%s %f','whitespace',' ');
new_line24             = fgets(fid2);
[dummy24,yllcorner2]   = strread(new_line24,'%s %f','whitespace',' ');
new_line25             = fgets(fid2);
[dummy25,cellsize2]    = strread(new_line25,'%s %f','whitespace',' ');
new_line26             = fgets(fid2);
[dummy26,NODATA_value] = strread(new_line26,'%s %f','whitespace',' ');

for i=1:1:nrows2
  new_line             = fgets(fid2);
  elev2(i,:)           = strread(new_line);
end

% Now generate Easting, Northing vectors:

for i=1:1:ncols2
  easting2(i)           = xllcorner2+(i-1)*cellsize2;
end

for i=1:1:nrows2
  northing2(i)          = yllcorner2+(nrows2-i)*cellsize2;
%  northing2(i)          = yllcorner2+(nrows2-(i-1))*cellsize2;
end

% Round off lidar easting, northing values:

easting2  = round(easting2);
northing2 = round(northing2);

% Radar DEM next:

fid1                  = ...
     fopen('./DEM_radar_200m_wgs84_utm57s.dat','r');
 
null_line              = '';
breaktest              = 0;
i                      = 0;

new_line11             = fgets(fid1);
[dummy11,ncols1]       = strread(new_line11,'%s %f','whitespace',' ');
new_line12             = fgets(fid1);
[dummy12,nrows1]       = strread(new_line12,'%s %f','whitespace',' ');
new_line13             = fgets(fid1);
[dummy13,xllcorner1]   = strread(new_line13,'%s %f','whitespace',' ');
new_line14             = fgets(fid1);
[dummy14,yllcorner1]   = strread(new_line14,'%s %f','whitespace',' ');
new_line15             = fgets(fid1);
[dummy15,cellsize1]    = strread(new_line15,'%s %f','whitespace',' ');
new_line16             = fgets(fid1);
[dummy16,NODATA_value] = strread(new_line16,'%s %f','whitespace',' ');

for i=1:1:nrows1
  new_line             = fgets(fid1);
  elev1(i,:)           = strread(new_line);
end

% Now generate Easting, Northing vectors:

for i=1:1:ncols1
  easting1(i)           = xllcorner1+(i-1)*cellsize1;
end

for i=1:1:nrows1
  northing1(i)          = yllcorner1+(nrows1-i)*cellsize1;
%  northing1(i)          = yllcorner1+(nrows1-(i-1))*cellsize1;
end

% Find areas common to radar, lidar images:

for i=1:1:ncols1
  if easting1(i)<easting2(1)
    lo_cols1 = i+1;
  elseif easting1(i)<easting2(ncols2)
    hi_cols1 = i;
  end
end

for i=1:1:nrows1
  if northing1(i)>northing2(1)
    lo_rows1 = i+1;
  elseif northing1(i)>northing2(ncols2)
    hi_rows1 = i;
  end
end

j                   = 0;

for i=lo_cols:1:hi_cols1
  j                 = j+1;
  easting_lidar(j)  = easting1(i);
end

ncols               = j;
j                   = 0;

for i=lo_rows1:1:hi_rows1
  j                 = j+1;
  northing_lidar(j) = northing1(i);
end

nrows               = j;

% Create radar, lidar, ice thickness arrays for common areas:

for i=1:1:nrows
  for j=1:1:ncols
    elev_radar(i,j) = interp2(easting1,northing1,elev1, ...
                              easting_lidar(j),northing_lidar(i),'linear');
    elev_lidar(i,j) = interp2(easting2,northing2,elev2, ...
                              easting_lidar(j),northing_lidar(i),'linear');
  end
end

% Calculate RMS difference between lidar, radar elevations:

rms_difference = sqrt(mean(mean((elev_radar-elev_lidar).^2)))

% Now follow TG_DRIVING_STRESS to calculate surface slopes along central
% flowline.  First we'll open up the flowline coordinate file:

fid3      = fopen('./tg_flowline_interp_coords.dat','r');
                                                                                
breaktest = 0;
i         = 0;
                                                                                
while breaktest==0
  new_line          = fgets(fid3);
  if new_line~=-1
    i               = i+1;
    [east_flowline(i),north_flowline(i),dist_flowline(i)] = ...
      strread(new_line,'%f %f %f','delimiter',',');
    length_flowline = i;
  else
    breaktest       = 1;
  end
end

% Now pick out the portion of the flowline within the area common to the 
% radar and lidar images:

j = 0;

for i=1:1:length_flowline
  if ((east_flowline(i)>=easting_lidar(1))& ... 
      (east_flowline(i)<=easting_lidar(ncols))& ... 
      (north_flowline(i)<=northing_lidar(1))& ... 
      (north_flowline(i)>=northing_lidar(nrows)))
    j = j+1;
    east_flow_comm(j)  = east_flowline(i);
    north_flow_comm(j) = north_flowline(i);
    dist_flow_comm(j)  = dist_flowline(i);
    n_flow_comm        = j;
  end
end

% Calculate easting, northing components of radar and lidar surface 
% slopes.  We'll be using an averaging scheme as follows:
%
% Slope = ((H(i+2)+H(i+1))/2)-(H(i-1)+H(i-2))/2)/(3*cellsize)
                                                                                
radar_slope_east         = zeros(nrows,ncols);
radar_slope_north        = zeros(nrows,ncols);
lidar_slope_east         = zeros(nrows,ncols);
lidar_slope_north        = zeros(nrows,ncols);
                                                                                
for i=3:1:nrows-2
  for j=3:1:ncols-2
    radar_slope_east(i,j)  = ...
      ((radar_elev(i,j-1)+radar_elev(i,j-2))/2 ...
      - (radar_elev(i,j+2)+radar_elev(i,j+1))/2)/(3*cellsize1);
    radar_slope_north(i,j) = ...
      ((radar_elev(i+2,j)+radar_elev(i+1,j))/2 ...
      - (radar_elev(i-1,j)+radar_elev(i-2,j))/2)/(3*cellsize1);
    lidar_slope_east(i,j)  = ...
      ((lidar_elev(i,j-1)+lidar_elev(i,j-2))/2 ...
      - (lidar_elev(i,j+2)+lidar_elev(i,j+1))/2)/(3*cellsize1);
    lidar_slope_north(i,j) = ...
      ((lidar_elev(i+2,j)+lidar_elev(i+1,j))/2 ...
      - (lidar_elev(i-1,j)+lidar_elev(i-2,j))/2)/(3*cellsize1);
  end
end
                                                                                
radar_slope_mag = sqrt(radar_slope_east.^2+radar_slope_north.^2);
lidar_slope_mag = sqrt(lidar_slope_east.^2+lidar_slope_north.^2);

radar_slope_flowline = interp2(easting_lidar,northing_lidar,elev_radar, ...
                       east_flow_comm,north_flow_comm,'linear');
lidar_slope_flowline = interp2(easting_lidar,northing_lidar,elev_lidar, ...
                       east_flow_comm,north_flow_comm,'linear');

%figure
%hold
%plot(dist_flow_comm,radar_slope_flowline,'ro-')
%plot(dist_flow_comm,lidar_slope_flowline,'bo-')
%xlabel('Along-Flowline Distance, km')
%ylabel('Surface Slope; r=radar; b=lidar')
%title('Radar and Lidar Surface Slope Comparison, 600m Smoothing')

% Now write results to a file:

fid11 = fopen('tg_lidar_comp.dat','w');

for i=1:1:n_flow_comm
  fprintf(fid11,'%d %d %d\n',dist_flow_comm(i),radar_slope_flowline(i), ...
                             lidar_slope_flowline(i));
end
