Feb 8, 2008

Import High Resolution DEM (2)

We are talking about more general case here: DEM is processed in ArcGIS (e.g. Mosaic, Extract by Mask, etc.). And it is ready for computation in Mathematica. The DEM can be exported as Float32 format, imported with the procedure in the previous post.
In this post, the DEM is exported as ARCGrid ASCII format, the whole DEM is in one text file.

data = Import["DEM.txt", {"Lines", All}];

(* the first lines are metadata *)

meta = data[[1 ;; 6]]

{"ncols 1703", "nrows 1200", "xllcorner
-86.528233463385", "yllcorner 41.430023494444", "cellsize
0.00027777777779994", "NODATA_value -9999"}

(* drop the metadata *)
data = Drop[data, 6];

(* convert string to data *)
data = Map[ToExpression, StringSplit[#]] & /@ data;

The next step is that we will make the output of ReliefPlot match the other geographic data.

1. Deal with NOData: PlotRange->{zMin,zMax}

2. Put the image into right geographic zone: DataRange-> {{Lon_Min, Lon_Max}, {Lat_Min, Lat_Max}}

ReliefPlot[data, DataReversed->True, DataRange-> {{-86.5282, -86.0552}, {41.4300, 41.7634}}, PlotRange -> {199.0, 283.0}]

Then tested it with the map created in another post.

Note: ReliefPlot needs an option: NoData->{}


Pooven said...

An interesting post, and quite useful! I'm trying to visualize the movement of air over a DEM using Mathematica. I understand that I'd need to use VecterFieldPlot3D but I'm uncertain as to how I'd represent the DEM in the 3D space? Do you have any ideas? I'm using the Show[] command which disallows the combining of my 3D vector space to the seemingly 2D DEM space.

I know that a DEM is meant to represent height given a coordinate, however when I open the DEM in Mathematica, the height is represented by colour instead of a 3D mesh. I'd probably need to feed the DEM in a different format? I'm uncertain as to the best way to do any of this. I have my ideas but I don't want to create something complex when the solution is one line :)

Anonymous said...

Thanks! Just stumbled across this post which has helped a lot. I needed to manipulate grid values from arc using mathematica functions which were too cumbersome to implement via map algebra.

In order to get the new grid values back to arc I copied the original header file and the Export command

Export[outFileFlt, FltNewValues, {"Binary", "Real32"}];

which worked fine. The only thing about the Float format is that arc geo-referencing info is lost. Do you have any tips for Importing / Exporting ESRI Raster formats directly in MMA?