1
$\begingroup$

I want to use .DEM files with GeoGraphics[]. For example, adding labels, or changing projections.

In my attempt, I import the example data from Wolfram.

Import["http://exampledata.wolfram.com/hailey.dem.gz", "ColorFunction" -> "DarkTerrain"]
data = Import["http://exampledata.wolfram.com/hailey.dem.gz", "Data"]
projection = Import["http://exampledata.wolfram.com/hailey.dem.gz", "Projection"]
spatialrange = Import["http://exampledata.wolfram.com/hailey.dem.gz", "SpatialRange"]

enter image description here

Adapted from this question, I generate the spacing of each pixel,

{rows, cols} = Dimensions@data
{rd, cd} = Differences /@ spatialrange
{rd, cd}/{rows, cols}
rds = spatialrange[[1, 1]] + Flatten@Accumulate[Table[rd/rows, rows]] + rd[[1]]/2;
cds = spatialrange[[2, 1]] + Flatten@Accumulate[Table[cd/cols, cols]];

I now try to convert this to geodata.

geodata = MapIndexed[GeoGridPosition[{cds[[#2[[2]]]], rds[[#2[[1]]]]}, projection] -> #1 &, data, {2}] // Flatten;

Which fails with invalid projection specifications:

GeoGridPosition::invproj: Invalid projection specification {None, "ReferenceModel" -> {Quantity[6.3565838*^6, "Meters"], Quantity[6.3782064*^6, "Meters"]}}.

Is there a smarter way to do this? Why is the projection not valid?

$\endgroup$

1 Answer 1

0
$\begingroup$

A solution is

geodata=MapIndexed[GeoGridPosition[{QuantityMagnitude[cds[[#2[[2]]]]],QuantityMagnitude[rds[[#2[[1]]]]]},{"Mercator","ReferenceModel"->{6378137,6378137}}]->#1&,data,{2}]//Flatten;

Which can be combined with GeoGraphics

GeoRegionValuePlot[MapAt[GeoPosition,geodata[[;;;;100]],{All,1}],GeoBackground->None,ColorFunction->GrayLevel]

and

enter image description here

It is very slow and not elegant. I might of been better creating some graphic that is stretched etc.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.