9
$\begingroup$

enter image description here

enter image description here

As a test, I'm trying to create a signed coastal distance map of the world (assumed perfectly spherical), using the first 1800x900 image above, which I converted to the 2nd image above using MorphologicalPerimeter, and was wondering if there was a more efficient approach, since I plan to use larger images to do something similar later.

My full code is below (also at https://github.com/barrycarter/bcapps/tree/master/MAPS/coastal.m), but my approach is to project coastal boundary points to 3D, create a region from them, and use the fast RegionDistance function.

Is there a better/faster way to do this? I know Mathematica has lots of Geo functions, but they seem to run slower and areharder (for me) to use unless I'm missing something.

(* Open the image, and create a boundary image *) image = Import["coastal.png"]; boundary = MorphologicalPerimeter[image]; (* use the image size to determine latitude/longitude of each pixel *) {height, width} = Take[Dimensions[ImageData[image]], 2] lats = N[Table[90-180*(i-1/2)/height, {i, 1, height}]]; lngs = N[Table[360*(i-1/2)/width-180, {i, 1, width}]]; (* perhaps excessive, but precomputing sins and cosines for efficiency *) colats = Map[Cos, lats*Degree]; colngs = Map[Cos, lngs*Degree]; silats = Map[Sin, lats*Degree]; silngs = Map[Sin, lngs*Degree]; (* convert latitudes and longitudes to unit sphere *) pts = Table[{colats[[i]] * colngs[[j]], colats[[i]]*silngs[[j]], silats[[i]]}, {i, 1, height}, {j, 1, width}]; (* select lit boundary pixels, create 3D region function *) litPixels = Position[ImageData[boundary], 1]; litPixels3D = Table[pts[[i[[1]],i[[2]]]], {i, litPixels}]; regionDistance = RegionDistance[Point[litPixels3D]]; (* apply regiondistance function to all points *) distances = Map[regionDistance, pts]; (* convert linear distances to geodesic distances in km *) lin2geo[d_] = 6371.009*2*ArcSin[d/2]; distancesGeo = Map[lin2geo, distances]; (* sign distances by determining if pixel was "lit" in original image *) imageData = ImageData[image]; distancesGeoSigned = Table[ distancesGeo[[i,j]]*If[imageData[[i,j]] == 1, -1, 1], {i, 1, height}, {j, 1, width}];
Improve this question
$\endgroup$

1 Answer 1

Reset to default
6
$\begingroup$

The DistanceFunction option in Nearest works with GeoDistance, effectively computing in 3D on the Earth's surface for us.

So we can

  • convert each pixel to a GeoPosition
  • identify the boundary GeoPosition objects
  • compute the unsigned geo distance field using Nearest
  • use the filled input image to infer the sign

Input data:

image = Import["https://i.sstatic.net/wAhoi.png"]; boundary = MorphologicalPerimeter[image]; {w, h} = ImageDimensions[image];

Nearest function:

toGeoPos = { Rescale[#2, {1, h}, {-90.0 + 0.5/h, 90.0 - 0.5/h}], Rescale[#1, {1, w}, {-180.0 + 0.5/w, 180.0 - 0.5/w}] }&; wpx = PixelValuePositions[boundary, 1]; nf = Nearest[Thread[GeoPosition[Transpose[toGeoPos @@ Transpose[wpx]]]] -> "Distance"];

Unsigned distance field (geo distance metric space):

LaunchKernels[]; dt = Developer`ToPackedArray[ParallelTable[ Developer`ToPackedArray[ QuantityMagnitude[ nf[GeoPosition[Thread[toGeoPos[Range[w], y]]]][[All, 1]], "Kilometers" ], Real ], {y, h, 1, -1}, Method -> "FinestGrained" ], Real]; // AbsoluteTiming 
{42.8758, Null}

Signed distance field:

sgn = ImageData[1 - 2image]; sdt = sgn * dt;

Colorized and legended image:

colorized = ImageMultiply[ Colorize[ImageAdjust[Image[sdt]]], ColorNegate[boundary] ]; legend = BarLegend[ {Blend["DarkRainbow", Rescale[#, MinMax[sdt]]]&, MinMax[sdt]}, LegendLabel -> Style["Signed Distance (km) ", GrayLevel[0.2], 10] ]; Legended[ Image[colorized, ImageSize -> 0.5w], Magnify[legend, 1.25] ]

enter image description here

Or stylize the ocean and land differently (harder to interpret the data):

landdf = ImageMultiply[ Colorize[ImageAdjust[ImageAdjust[Image[Ramp[-sdt]]], {0,0,.5}], ColorFunction -> "AlpineColors"], image ]; oceandf = ImageMultiply[ Colorize[ImageAdjust[ImageAdjust[Image[Ramp[sdt]]], {0,0,.65}], ColorFunction -> "DeepSeaColors"], ColorNegate[image] ]; ImageAdd[landdf, oceandf, boundary]

enter image description here

We can also look at a relief plot or contours:

autolegend = BarLegend[Automatic, LegendLabel -> Style["Signed Distance (km) ", GrayLevel[0.2], 14]]; ReliefPlot[sdt, DataReversed -> True, Frame -> False, PlotLegends -> autolegend]

enter image description here

ListContourPlot[Reverse[sdt], AspectRatio -> Automatic, Frame -> False, PlotLegends -> autolegend]

enter image description here

And we can look at these signed distance fields in 3D:

im = ImageAdd[landdf, oceandf, boundary]; Graphics3D[ {Texture[im, "Spherical"], Black, Sphere[{0, 0, 0}]}, Boxed -> False, ViewPoint -> {-2.12, -2.5, 0.83}, ViewVertical -> {0, 0, 1}, Lighting -> "Neutral" ]

Lastly, we should note our sdf is only as accurate as our input data. For example since we're missing some small islands, this field misidentifies point nemo.

Correct point nemo (in yellow below): GeoPosition[{-48.8767, -123.393}]

Our computed point (as the red × below):

GeoPosition[toGeoPos @@ ({0,h} + {1,-1}Reverse[Position[sdt, Max[sdt]][[1]]])] 
GeoPosition[{-25.1278, -108.56}]

However we're only ~200 miles off from the Eurasian pole of inaccessibility. Wikipedia's ground truth in green and our computed value in purple:

enter image description here

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.