% gunzip the South San Francisco DEM file to a temporary directory
filenames = gunzip('sl3d_sanfranciscos.dem.gz', tempdir);
demFilename = filenames{1};

% read the 1:24,000 DEM file
% standardize missing values to NaN
deminfo = georasterinfo(demFilename);
Z = standardizeMissing(readgeoraster(demFilename), deminfo.MissingDataIndicator);

% delete the temporary gunzipped file
delete(demFilename);

demdata=Z;
[xdim, zdim] = size(demdata);
xspace = 30; % scaling in meters for x dimension
zspace = 30; % scaling in meters for z dimension
% reshape the data into a one-dimensional array
demdata = demdata(:);

% bring up template
myworld = vrworld('vr_template_terrain.x3d');
% open the virtual world
open(myworld);
% create a handle to a node VRTerrain, the node that will contain the DEM data
Terrain_node = vrnode(myworld,'VRTerrain');

% create a child of VRTerrain - shape
newShape = vrnode(Terrain_node, 'children', 'Terrain_Shape', 'Shape');
% create appearance field for the shape
newAppear = vrnode(newShape, 'appearance', 'Terrain_Appearance', 'Appearance');
% create material field for the appearance
newMat = vrnode(newAppear, 'material', 'Terrain_Material','Material');
% assign properties for the material field
newMat.ambientIntensity = 0.25;
newMat.diffuseColor = [0.9 0.6 0.6];
newMat.shininess = 0.078125;
newMat.specularColor = [0.0955906 0.0955906 0.0955906];
% create geometry field for the shape
newEGrid = vrnode(newShape, 'geometry', 'DEM_EGrid','ElevationGrid');
% assign properties for the geometry field - use DEM data
newEGrid.creaseAngle = 3.14;
newEGrid.xDimension = xdim;
newEGrid.zDimension = zdim;
newEGrid.xSpacing = xspace;
newEGrid.zSpacing = zspace;
newEGrid.height = demdata;
newEGrid.ccw = 'TRUE';
% This setting will make the terrain surface visible from both sides
newEGrid.solid = 'FALSE';

% terrain elevation is used to color the image
cmap = demcmap(Z, 256);

% create texture subdirectory of the current directory
% output arguments used only to avoid warning message when the directory
% already exists
[~, ~] = mkdir('texture');

% scale the height values to use the full colormap range
% scaling relies on the fact that this terrain begins at zero height
Zscaled = Z .* (size(cmap,1)-1) ./ max(Z(:));

% save the texture into PNG image in the texture subdirectory
% rotate the image left to match image orientation needed in X3D model
% elements of Zscaled represent indices into cmap
imwrite(rot90(Zscaled), cmap, 'texture/sanfrancisco_elev.png');

newTexture = vrnode(newAppear, 'texture', 'Terrain_texture','ImageTexture');
newTexture.url = 'texture/sanfrancisco_elev.png';

% copy b747.x3d from /sl3ddemos to the current directory
pt = fileparts(which('vrterrain_simple.m'));
copyfile(fullfile(pt, 'b747.x3d'), pwd, 'f');

% create a new Transform node, called "Boeing"
plane = vrnode(myworld, 'Boeing', 'Transform');
plane_inline = vrnode(plane, 'children', 'Boeing_Inline', 'Inline');
% a simple model of Boeing is prepared in the /sl3ddemos directory
plane_inline.url='b747.x3d';

ypeak = max(Z(:));
[xmax, zmax] = find(Z==ypeak);
% use the first peak, if more vertices have the same maximum height
% convert matrix indices to meters in x and z directions in X3D
xpeak=xspace*(xmax(1)-1);
zpeak=zspace*(zmax(1)-1);

plane.translation = [xpeak ypeak+200 zpeak];
% scale the size of the airplane by a factor of 20, so that it
% is visible in the virtual scene without any extra zooming
plane.scale = [20 20 20];

%
% The triad consists of 3 lines (1 meter long) running from one vertex
% along the x, y and z directions. The lines are coloured as follows:
% +x - red
% +y - green
% +z - blue
%
% If the triad is included in the scene at the top level of the scene
% hierarchy, it denotes the global scene coordinates.
% If it is included as a child of a transform node, it denotes the local
% coordinate system (orientation) of that node in the scene.

% add the triad to the scene at the top level of the hierarchy
triad = vrnode(myworld, 'Triad1', 'Triad');
% scale the size of the triad so that it is visible
% in the virtual scene without any extra zooming
triad.scale = [xdim*xspace/8 min(xdim*xspace/8, zdim*zspace/8) zdim*zspace/8];
% position the triad at the center of Boeing 747
triad.position=[xpeak ypeak+200 zpeak];

myworld.World_Info.title = 'B747 Flying over the San Francisco Area';

save(myworld, 'vrterrain_sanfrancisco.x3d');

close(myworld);
delete(myworld);

% This is how to open a X3D file in the external viewer:
% vrview('terrain.x3d', '-web');

% This is how to open a X3D file in the internal viewer:
% vrview('terrain.x3d', '-internal');

% This is how to open the X3D file in the default viewer:
createdworld = vrview('vrterrain_sanfrancisco.x3d');
% Set Antialiasing on to smooth the terrain texture
myfig = get(createdworld,'Figures');
set(myfig, 'Antialiasing', 'on');

% clear all used variables
clear Terrain_node Z Zscaled cmap createdworld demFilename demdata ex filenames ...
  id lat lon mess myfig myworld newAppear newEGrid ...
  newMat newShape newTexture nm ok plane plane_inline pt triad ...
  ve xdim xmax xpeak xspace ypeak zdim zmax zpeak zspace