function [cn, X, Y, data] = load_dem(file_name)

fid=fopen(file_name);
    firstline = fgetl(fid);
    ncols = str2num(firstline(7:length(firstline)))
    secondline = fgetl(fid);
    nrows = str2num(secondline(7:length(secondline)))
    thirddline = fgetl(fid);
    xllcorner = str2num(thirddline(11:length(thirddline)))
    fourthline = fgetl(fid);
    yllcorner = str2num(fourthline(11:length(fourthline)))
    fifthline = fgetl(fid);
    cellsize = str2num(fifthline(10:length(fifthline)))
    sixthline = fgetl(fid);
    NODATA_value = str2num(sixthline(14:length(sixthline)))
    data = fscanf(fid, '%g', [ncols nrows]);
    size(data)
fclose(fid);

cn = length(find(data==NODATA_value)); %count nulls
data(find(data==NODATA_value)) = NaN;
x = xllcorner:cellsize:(xllcorner+(ncols-1).*cellsize);
yy = (yllcorner+(nrows-1).*cellsize):-cellsize:yllcorner;
[X, Y]=meshgrid(x,yy);


figure(1)
clf
plot(xllcorner, yllcorner, 'k.')
hold on
imagesc(x,yy,data')
colormap jet
colorbar
axis equal
xlabel('distance easting (m)')
ylabel('distance northing (m)')
%plot(X,Y,'k.')



figure(2)
clf
hold on
imagesc(x,yy,gradient(data'))
hold on
colormap copper
axis equal
xlabel('distance easting (m)')
ylabel('distance northing (m)')