Importing the DELA-4099.xyz terrain data

In this worksheet, we model importing the DELA-4099.xyz terrain data into Mathematica and using it to construct a function which gives terrain height (in meters) as a function of latitude and longitude. I got this data from the NOAA website "Design-A-Grid" custom grid generator at 

http://www.ngdc.noaa.gov/mgg/gdas/gd_designagrid.html#

To get it, I entered the lat and long of the corners of the grid I was interested in (110 W and 120 W, 40 N and 30 N, obtained from browsing around on Google Earth), selected XYZ data,  no header, comma, and unchecked "Omit Empty Grid Cells". I got a .tgz file with a bunch of descriptive stuff and the file DELA-4099.xyz in it. Looking at the number of grid cells, I saw that my data would be 601 x 601.

To get this into Mathematica, we use Mathematica's Import command to import the data into a Mathematica table.

?Import

We'll have to find the file DELA-4099.xyz on our system. Here's a list of directories where Mathematica will look for it:

$Path

Take a moment and move the file DELA-4099.xyz to one of these locations. Now we can import the file. Since we selected "Comma" on the NOAA page, our data is comma-separated, so we give "CSV" (Comma Separated Values) as the format. Note that there is an "XYZ" format, but Mathematica won't parse the file generated by the NOAA page directly as XYZ.

RawHeightData = Import["DELA-4099.xyz","CSV"];

We want to massage this data so that we get heights on an evenly spaced grid in the format

{{z11,z12,z13,z14... },{z21,z22,z23,z24,... },                }

where each sublist has the same x coordinate and the y coordinates vary across the sublist (weird, but that's the way it works).
To do that, we want to collect all of the data points which have the same x value and make them sublists. We'll start by sorting the data:

SortedHeightData = Sort[RawHeightData];

Now we can generate the list of lists (knowing that we have a 601 x 601 dataset) using

GridData = Partition[SortedHeightData,601];

We now generate the table by taking the 3rd part of each sublist

HeightGrid = GridData[[All,All,3]];

And finally we can come up with an InterpolatingFunction for the height data using ListInterpolate:

z = ListInterpolation[HeightGrid,{{-120,-110},{30,40}}]

We get a perfectly good Mathematica function z, whose inputs are in lat and long and whose output is elevation in meters. Let's plot that a couple of ways. First, here's the Rockies as a 3D plot. Why are the latitiudes negative? Because they are measured west of 0 degrees (Greenwich) rather than east. So the -120 latitude is the most westerly line on the map while the -110 latitude line is the farthest east.

Plot3D[z[x,y],{x,-120,-110},{y,30,40},ColorFunction->ColorData["LightTerrain"]]

Next, here's the contour plot:

ContourPlot[z[x,y],{x,-120,-110},{y,30,40},ColorFunction->ColorData["LightTerrain"]]