# How can I extract the height of a DEM point?

Posted 9 years ago
13924 Views
|
11 Replies
|
4 Total Likes
|
 I have a DEM that makes a fine picture when I use the Input function Import[ "C:\\GIS \ Louisiana\\30091472_Plaquemines_NE\\3009147ne_dem\\dem_3009147ne.dem"] I wish to test/manipulate the elevations at a set of sample locations. I've looked at the elements imported but cannot figure out how to query the elevation data for one location (pixel). Import[ "C:\\GIS \ Louisiana\\30091472_Plaquemines_NE\\3009147ne_dem\\dem_3009147ne.dem","Elements"] Any guidance very appreciated.JAC
11 Replies
Sort By:
Posted 9 years ago
 I imported that DEM file in R using the following commands: # DEM file # Set working directory setwd("C:\\users\\jbaldwin\\Downloads") # Load library library(rgdal) # Input data g = readGDAL("dem_3009147ne.dem") # Get coordinates xy = coordinates(g) # Show the structure of g str(g) So getting the coordinates (in UTM's) can be done by running R from Mathematica. Alternatively by comparing what R gets and the elements that Mathematica imports, one can determine the UTM's of any element of the data matrix with (* Get necessary elements of DEM file *) elev = Import["C:\\Users\\jbaldwin\\Downloads\\dem_3009147ne.dem", "Data"]; dim = Import["C:\\Users\\jbaldwin\\Downloads\\dem_3009147ne.dem", "Dimensions"] resolution = Import["C:\\Users\\jbaldwin\\Downloads\\dem_3009147ne.dem", "SpatialResolution"] boundingBox = Import["C:\\Users\\jbaldwin\\Downloads\\dem_3009147ne.dem", "SpatialRange"] followed by (* Get elevation (in meters) and coordinates (in UTMs) for a \ particular element *) i = 1291 j = 350 elev[[i, j]] utmx = boundingBox[[2, 1]] + i resolution[[1]] utmy = boundingBox[[1, 1]] + j resolution[[2]] or a something like a Table that gets all of the UTM's for each row and each column. You'll be getting UTM's rather than latitudes and longitudes. Also, the above R code ends up with the elevation in feet. Either there's some flag in the file that says it's in feet and R keeps it in feet or Mathematica always converts to meters for this particular file format irrespective of that flag or certainly some other rule. I don't know.
Posted 9 years ago
 Thank you, Dr. Purcell. I am fairly confident of calculating the Lat, Lon as needed, but was struggling to get the elevation. I will practice with your contribution to this neophyte. Your suggestions look promising. I haven't yet seen the Mathematica documentation for "DEM", which disturbs me because of all the help file reading I've been doing.Dr. Reisse, I will have to re-read the DEM specifications to be sure, but I believe the file format locates one corner and defines the spacing of the points in Lat, Lon and projects it in UTM. One result is a thin wedge on each edge that contains NULL values, because the projection is planer and the point spacing is geodetic.Dr. Meyer of U.Conn. has given me advice outside of this thread for which I am grateful. His advice has developed into a bit of a collaboration (his programming skill seems a couple of orders of magnitude above my own). Perhaps we will have something to share in the near future.Thank you all for your contributions!JAC
Posted 9 years ago
 Interpolating the data set contained in the DEM file is quite straightforward (see the earlier post above) as is getting the other data elements contained in the file. The key questions that are being raised here are first, how to get the actual latitude and longitude range from the file. The data elements do not appear to contain this. The second issue is, since the data represents a grid for a particular map projection, how to invert that to actually get the latitude the height for specific latitude and longitude points once those ranges are known.
Posted 9 years ago
 If you search in the Mathematica documentation for "DEM" you will see options for reading the data and metadata, which you will probably know how to interpret. In this example I am just using pixel count for the coordinates of the x and y axes, since I don't know how to convert to lat and long. Tooltip is used to display the elevation, when the cursor is put over a contour line. This is not quite what you want, but close. file = "/Users/christopherpurcell/Desktop/dem_3009147ne.dem"; data = Import[file, "CoordinateSystemInformation"] bounds = Import[file, "SpatialRange"] (* should be {lat, long} in dec degrees but must \ be something else, maybe meters *) data = Import[file, "Data"]; size = Import[file, "Dimensions"] (* data array size *) interp = ListInterpolation[Transpose[data]]; ContourPlot[ Tooltip[interp[x, y]], {x, 1, size[[2]]}, {y, 1, size[[1]]}, Contours -> Range[10], PlotPoints -> 50] The next step might be to put a ReliefPlot of the data in a Manipulate, and return the elevation value that corresponds to the cursor location. Something along these lines, where if you click on a point, you get a display of the elevation in meters. (I think). Manipulate[ ReliefPlot[ data , FrameTicks -> True, PlotLabel -> ToString[Join[x, {interp[Apply[Sequence, x]]}]], ColorFunction -> "LightTerrain", PlotLegends -> Automatic], {{x, {1, 1}}, Locator}] Next you would figure out the lat and long and use DataRange option to get the axes scaled correctly.
Posted 9 years ago
 I am trying to see if maybe I can kludge a solution. The following demonstrates the ability to outline a parish (county) as a polygon. GeoGraphics[{GeoStyling["StreetMap"], Polygon[Entity["AdministrativeDivision", {"EastBatonRouge", "Louisiana", "UnitedStates"}]]}, GeoBackground -> None] If one double clicks and double clicks again in the center of the image all the nodes of the polygon are highlighted. What I haven't been able to do is to tease out the?coordinates of the nodes. With a matrix of the coordinates I may be able to recursively test inside the polygon for the characteristics I am looking for. Then I can repeat for other parishes and counties.I've played with Entities and GeoEntities without any luck, so far.My earlier hope was to simply use a matrix of the height values and coordinates (bounded by the polygon) for testing. If one knows a coordinate the height is available. Taking many tthousands ?of points is beyond manual input methods.JAC
Posted 9 years ago
 Thanks for the shout-out, Dr. Reiss.I'm admittedly new to Mathematica and when I saw easy to use functions like Import["...*.dem"] and GeoGraphics and GeoStyling and DensityPlot, that obviously access the height information to choose colors, I thought I would be able to interrogate the same values that give colors to give me a numerical value, but I've not been able.JAC
Posted 9 years ago
 Do any of the WRI developers on the forum have any advice for Mr. Cavell? The main issue is, given a DEM file, read it in and generate a function height[latitude,longitude] based on an interpolation of the supplied grid data in the DEM file. The two things that need to be extracted and analyzed for this area) the latitude and longitude range of the DEM array data b) any transformations that need to be done on the data based on the map projection indicated in the DEM fileIt would be a good addition to DEM import if such a function could be generated as one of the import Elements for DEM
Posted 9 years ago
 It seems I've hit a dead-end. I just wish I was able to take apart the pieces that are obviously present to enable the graphical representation, but it seems "over my head." Perhaps reverting to something other than Mathematica will provide an answer.JAC
Posted 9 years ago
 I wish to compliment Dr. Reiss for his help. I contacted him off-line because I had a follow-up question that involved a brand name other than Wolfram and didn't want to be impolite to our host. He has suggested a few ideas that seem to be getting me closer to a possible solution. file = "...local folders.../dem_3009147ne.dem" Import[file, "ReliefImage"] interp = ListInterpolation[Transpose@data] interp["Domain"] {{1., 1225.}, {1., 1405.}} DensityPlot[interp[x, y], {x, 1, 1225}, {y, 1, 1405}, PlotPoints -> 100] He then suggested I revisit the GeoElevationData function, which I did. Providing a lat, log it returns a value for height, but I am wary of it because I haven't found a quotable statement of the source of that elevation result. One reason I was investigating the DEM was I know its pedigree and can stand behind the results.A point on the peninsula looking crook in the river is 30°19'48.13"N lat, 91° 9'38.35"W lon (30.33003611, -91.16065278) GeoElevationData[{30.33003611, -91.16065278}] returns: Quantity[55.7743, "Feet"] which may be an OK value but I can't be sure. The grid on which it is based is pretty coarse compared to the DEM and when I checked specific sites in urban areas the discrepancies cause some concern.JAC
Posted 9 years ago
 Thanks for looking into this. The subject file (dem_3009147ne.dem) is attached. I have also tried to do something similar with embedded elevation models, also with no joy so I am trying with a local file. An example of a previous attempt is: GeoGraphics[{EdgeForm[Black],GeoStyling["ReliefMap"],Polygon[=East Baton Rouge Parish]}, GeoBackground -> None] (The = in front of "East Baton Rouge Parish" is ctrl-=.)JAC Attachments:
Posted 9 years ago
 Could you attach your example file so that we have a chance to play with it and see what we can suggest?