%%% This script reads in a dem from NASA SRTM data repository here: 
%%%    ftp://e0srp01u.ecs.nasa.gov/srtm/
%%% Files are located here: 
%%%    ftp://e0srp01u.ecs.nasa.gov/srtm/version2/SRTM3
%%% 
%%% Note: a volcano name must be provided along with the parent directory where
%%% where the volcano dem is located.  Cropping option selects a square
%%% region around the point of interest (often provided as the summit or
%%% vent of a volcano).
%%%
%%% lat_range and lon_range provides options for cropping an original 1
%%% degree by 1 degree grid.  Integer arguments corresponding to dem
%%% latitude and longitude spacing, which is three seconds for
%%% international maps.  Default provides +/-5 minutes of latitude and
%%% longitude
%%%
%%% Jeff Johnson - June 23rd 2009 - GEOP572



function [dem_zoom,lats_zoom,lons_zoom] =  mappa(volcano,lat_range,lon_range)
%mapsize = 12; % kilometers from map center

cd '/Users/andrea/Desktop/SEGY_tungu09' % directory of SRTM volcano maps
%cd(volcano)

if nargin<2;
    lon_range = 100; % +/-5 minutes longitude default
    lat_range = 100; % +/-5 minutes latitude default
end

switch volcano
    case 'TUNGURAHUA'
        %lat = -1.467; lon = -78.442; % Tungurahua summit
        lat = -1.43427; lon = -78.48844; % TLEF
        filenames = ['S02W079.hgt'];
        latrange = [-2 -1];
        lonrange = [-79 -78];
end

%%% Note that some volcanoes are located at the boundary of several maps
%%% such that multiple maps must be stitched together.
mapdata = [];
for i=1:size(filenames,1)
    filename = filenames(i,:)
    fid = fopen(filename,'r','ieee-be');
    mapdem = fread(fid,'int16');
    fclose(fid);
    npts = sqrt(length(mapdem));
    mapdem = flipud(reshape(mapdem,npts,npts)');
    mapdem(mapdem<0) = nan;
    if i==1 % only one map grid is red in
        mapdata = mapdem;
    elseif diff(filenames(:,7))==0 % case where SRTM maps are aligned vertically
        mapdata = [mapdata; mapdem(2:npts,:)];
    elseif diff(filenames(:,3))==0 % case where SRTM maps are aligned horizontally
        mapdata = [mapdata mapdem(:,2:npts)];
    end
end
dem = mapdata;
lats = linspace(latrange(1),latrange(2),size(mapdata,1));
lons = linspace(lonrange(1),lonrange(2),size(mapdata,2));

lat_index = find(min(abs(lats-lat))==abs(lats-lat));
lon_index = find(min(abs(lons-lon))==abs(lons-lon));
lats_zoom = lats(lat_index-lat_range:lat_index+lat_range);
lons_zoom = lons(lon_index-lon_range:lon_index+lon_range);
dem_zoom = dem(lat_index-lat_range:lat_index+lat_range,...
    lon_index-lon_range:lon_index+lon_range);
dem_zoom = interp2d(dem_zoom,lats_zoom,lons_zoom);