Converting BIL terrain data into an Octave/MATLAB array

This one was a real brain-buster. The USGS has an awesome app called EarthExplorer (here) that allows you to select a region (or regions), see what data is available for them, and then order the data. The application for ordering the data, BDA, gave me several issues when run on Windows 10. Switching to Ubuntu for using BDA completely fixed the problem.

Parsing the BIL files that were provided was a bit tricky in Octave. Plenty of resources are available for that kind of thing if you’re using Python, but I found myself creating an m-file from scratch to make this project work.

function  [lats longs alts ds_lats ds_longs ds_alts] = bil2mat()

% bil2mat.m
% Developed in Octave 4.4.1 by Dover Test Labs // doverslab.com
% 
% This function is not called with any parameters, rather it prompts the user to 
% select a BIL file. The BIL header file is assumed to be in the same folder.
% The BIL file is then parsed and output as 1xn vectors (lats,longs,alts).
%
% For display purposes, the vectors are downsampled and output as (ds_* vectors)
%


user_yes = 1;

% these are incremented versions of the lat/long/alt parameters
longi = [];
lati = [];
alti = [];

% these are the instantiations of the final output vectors
longs = [];
lats = [];
alts = [];

% these are the instantiations of the downsampled output vectors
ds_longs = [];
ds_lats = [];
ds_alts = [];

while user_yes == 1 % loop while the user continues to add additional BIL files
    [file path] = uigetfile

    hdr_str = strrep(file,'.bil','.hdr');

    [hdrs,;] = importdata([path hdr_str],'\t');
    
    num_hdrs = size(hdrs,1);
    
    for ii = 1:num_hdrs
      
      parm = strsplit(hdrs{ii,1});

   
      switch(parm{1})
        case "NCOLS"
          NCOLS = parm{2};
        case "NROWS"
          NROWS = parm{2};
        case "ULXMAP"
          ULXMAP = parm{2};
        case "ULYMAP"
          ULYMAP = parm{2};
        case "YDIM"
          YDIM = parm{2};
        case "XDIM"
          XDIM = parm{2};
        case "NODATA"
          NODATA = parm{2};
      endswitch
             
    end
    
    NROWS = str2num(NROWS)
    NCOLS = str2num(NCOLS)
    ULXMAP = str2num(ULXMAP)
    ULYMAP = str2num(ULYMAP)
    YDIM = str2num(YDIM)
    XDIM = str2num(XDIM)
    NODATA = str2num(NODATA)
      
    bfile = fopen([path file]);
    alt_vec = fread(bfile,inf,'int16');
    alt_mat = reshape(alt_vec,NROWS,NCOLS)'; 
    
    % Find indexes of NODATA points
    ND_inds = find(alt_mat == NODATA);
    Num_NDs = length(ND_inds);
    alt_mat(ND_inds) = 1500;
    
    ND_inds_v = find(alt_vec == NODATA);
    alt_vec(ND_inds_v) = 1500;    
    
    disp('Number of NoData Pts found:')
    disp(Num_NDs)
    

    % Calculate lats, longs, alts using interval from header file   
    longi = linspace(ULXMAP,ULXMAP+NCOLS*XDIM,NCOLS);
    lati = linspace(ULYMAP,ULYMAP-NROWS*YDIM,NROWS);
    alti = alt_vec;
    
    disp('size of alt_vec')
    disp(size(alt_vec))   
   
    lon_vec = repmat(longi,1,NCOLS);
    
    disp(size(lon_vec));
    
    lat_vec = repelem(lati,1,NROWS);
    
    disp(size(lat_vec));

    %%%%%% Downsampling for demo purposes %%%%%%%%
    ds_rate = 100;
    r_lin_iis = 1:ds_rate:NROWS;
    c_lin_iis = 1:ds_rate:NCOLS;
     
    long_ds = longi(r_lin_iis);
    lat_ds = lati(c_lin_iis);
    alt_ds = alt_mat(r_lin_iis,c_lin_iis);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    [X,Y] = meshgrid(long_ds,lat_ds);
    Z = griddata(long_ds,lat_ds,alt_ds,X,Y);
    hold on
    surf(X,Y,Z);

    btn = questdlg('Add more terrain?','Keep Going?');
    
    alts = [alts alt_vec'];
    longs = [longs lon_vec];
    lats = [lats lat_vec];
    
    ds_alts = alts(1:ds_rate:length(alts));
    ds_longs = longs(1:ds_rate:length(alts));
    ds_lats = lats(1:ds_rate:length(alts));
    
    %%% Scale graph based on Latitude
    lat_swath = max(ds_lats) - min(ds_lats);
    earth_rad = 6378137;
    circ = earth_rad*2*pi();
    lat_sw_m = circ*lat_swath/360;
    
    long_swath = abs(max(ds_longs) - min(ds_longs));
    long_circ = circ*cosd(mean(max(ds_lats) - min(ds_lats)));
    long_sw_m = long_circ*long_swath/360;
    
    alt_sw_m = max(max(ds_alts)) - min(min(ds_alts));
    
    zlim([0 0.25*mean([lat_sw_m long_sw_m])]);
            
    if strcmp(btn,'No')
      user_yes = 0;
    endif
 

 
 end
 
 hold off
    
endfunction