Different Axes Scales on 3D Graphics

Questions
,
1

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.

2

I have found a source of very accurate terrain elevation data for the Alpes, <snip>

Can you possibly share this source here in order to try a possible solution out prior to publication?

3

Hi, Udo

This is where I found the data…

Regards

Brad

4

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 content

ncols        3601
nrows        3601
xllcorner    5.999861111111
yllcorner    46.999861111111
cellsize     0.000277777778
NODATA_value -32768

but 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

<PAMDataset>
  <Metadata>
    <MDI key="AREA_OR_POINT">Point</MDI>
  </Metadata>
  <PAMRasterBand band="1">
    <UnitType>m</UnitType>
  </PAMRasterBand>
</PAMDataset>

testifying the unit is m, which should mean meters … something seems to elude me …

5

Hello, Udo

In 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?

Regards

Brad

6

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

7

Typing

vareyPlot["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 acount

bradVarey47N7E.png
bradVarey47N7E.png1000×1000 158 KB

first 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.

8

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}]