0
|
6624 Views
|
7 Replies
|
1 Total Likes
View groups...
Share

# Different Axes Scales on 3D Graphics

Posted 11 years ago
 I have found a source of very accurate terrain elevation data for the Alpes, and want to plot it using something like ListPlot3D.The data file is simply a list of 16bit signed integers which represent elevations in meters. There are approximately 3,000 rows of 3,000 columns, each row representing a given latitude, each column a given longitude. Points are 30 metres apart in each horizontal direction.I've got all the graphing working, but am now having a problem with scale. I want my X and Y axes to show lat and long in decimal degrees, and the Z axis to plot elevation in meters.But if you don't scale the axes, the view is distorted from how it would look naturally: a degree of longitude at 45 degrees latitude is much shorter in length than a degree of latitude; a meter does not scale directly to degrees.So I need to scale two of my axes as a function of the third. I know how to do the math, but I don't know how to tell, say, ListPlot3D, how to scale the axes. Of course I could convert degrees to meters, but then my axis labels will not show lat-long, which is what I want.Any help would be greatly appreciated.
7 Replies
Sort By:
Posted 10 years ago
 For a fast insight there is ReliefPlot Clear[elevationFile]  elevationFile = FileNameJoin[{NotebookDirectory[], "others", "elevation", "N47E006.hgt"}]; Clear[fileColumns, fileRows, arcSecondResolution] If[FileByteCount[elevationFile] < 3000000, {fileColumns = 1201;     fileRows = 1201; arcSecondResolution = 3}, {fileColumns = 3601;     fileRows = 3601; arcSecondResolution = 1}   ]; Clear[elevations]elevations = BinaryReadList[elevationFile, "Integer16", ByteOrdering -> 1];Clear[elevationsGrid]elevationsGrid = ArrayReshape[elevations, {fileRows, fileColumns}];ReliefPlot[elevationsGrid, AspectRatio -> 1./Cos[47. Degree], BoxRatios -> {Cos[47. Degree], 1, 10}]
Posted 10 years ago
 TypingvareyPlot["others/elevation/N47E006.hgt", 5.999861111111, 46.999861111111, 0.000277777778, 6371000.8/10(* 6371000.8 *), 20]gives that picture (mountains are a factor of 10 higher (or earth has 1/10 of its real radius)) and only every 20-th point is taken into acountfirst one displays the data on a sphere then one rotates to have the normal at 47N 6E pointing into {0,0,1} direction. Specifications in that area are for people who know it already. I did not figure out whether the lower-left corner point is the first data point in the file - as I presume - or not, and if not, where it is. In any case you might rework the notebook or discuss ways to fit it better to your needs. Attachments:
Posted 11 years ago
 Hello Brad, Does this help?Yes, a lot, thank you; if the week-end is a bit rainy I will have a look into it hopefully
Posted 11 years ago
 Hello, UdoIn terms of reading the file, I took a different way around, as I didn't try to convert it...elevationFile = "/Users/Brad/Documents/Wolfram/Misc/N45E006.hgt";setFileParameters;elevations =   BinaryReadList[elevationFile, "Integer16", ByteOrdering -> 1];(* First partition the whole file into rows and colums *)elevationsGrid = Partition[elevations, fileRows, fileColumns];The variables fileRows and fileColumns were set earlier by looking at the file size, as these files always have either 90 or 30 metre resolution for the same square degree. So a fragment of my function setFileParamters[] above does this...(* If the file is under 3MB, it is 3-arc-second resolution; otherwise \1-arc-second *)If[  FileByteCount < 3000000,   {   fileColumns = 1201;   fileRows = 1201;   arcSecondResolution = 3;   },   {   fileColumns = 3601;   fileRows = 3601;   arcSecondResolution = 1;   }  ];Note the "big-endian" order of the bytes, which is handled by the line above...BinaryReadList[elevationFile, "Integer16", ByteOrdering -> 1];Do all of that, and the array "elevations" (for a 1 arc-second file) contains values which are metre elevation points, derived from 16 bit integers, 3601 rows of 3601 columns, each corresponding to an arc-second point within the grid.You're correct about the altitude of Mont Blanc, but if that one highest point 4810 metres doesn't appear directly under one of the sample points, the maxium value of the whole array will be something slightly less, because the sampling would fall somewhere on the shoulder, and not the peak of the mountain. From memory, I think the highest point was 4806, or something close.Does this help?RegardsBrad
Posted 11 years ago
 Well, In[2]:= SetDirectory[ FileNameJoin[{NotebookDirectory[], "others", "elevation"}]] Out[2]= "N:\\Udo\\Abt_N\\others\\elevation"  In[3]:= elevD = Import["N47E006.hgt", "HGT"] During evaluation of In[3]:= Import::format: Cannot import data as HGT. >> Out[3]= $Failed In[8]:= Position[$ImportFormats, "HGT"]Out[8]= {}Mma does not read HGT as it stands now. So I used http://converter.mygeodata.eu/raster to convert into AAIGrid etc. pp. As I understand the file delivers a raster, X and Y extrema are given in the header while the height values are the contentncols        3601nrows        3601xllcorner    5.999861111111yllcorner    46.999861111111cellsize     0.000277777778NODATA_value -32768but then appear only three digit numbers, how are they separated to give heights which  need more than 3 digits in the length unit  (Mont Blanc is 4810 meters)? There is a XML file       Point        m  testifying the unit is m, which should mean meters ... something seems to elude me ...
Posted 11 years ago
 Hi, UdoThis is where I found the data...http://www.viewfinderpanoramas.org/dem3.htmlRegardsBrad
Posted 11 years ago
 I have found a source of very accurate terrain elevation data for the Alpes, Can you possibly share this source here in order to try a possible solution out prior to publication?