# Project map to a particular shape

I want to project a map to a particular shape, like a circle or ellipse. For example, let's imagine we want to project a map of the continental United states so that its outline is that of a circle. The basic strategy I have is illustrated below:

To project the outline is easy, the geographic center is at the center of the circle and the perimeter is divided into small sectors, perhaps 1-degree each. Then the points of the perimeter are equally divided within their sector and mapped to the corresponding sector on the perimeter of the circle.

To project the interior points of the map is the tricky part on which I need advice. In the image an example of such a point is shown as a blue dot. Then the procedure I have in mind for projecting it is as follows.

(1) Find the point on the geodesic perimeter of the United States which is closest to the point to be projected. This is indicated as a red dot in the illustrative image above.

(2) In the geodesic data (a sphere), draw an imaginary line to the red dot, the closest point.

(3) In the geodesic data (a sphere), draw two lines going outwards at 120 degrees in either direction relative to the red line. These two lines are shown in green. Find the points where these lines intersect the geodesic perimeter of the United States. These two points are shown as green dots.

(4) To project the point of interest, the blue dot, start at the red dot and move inwards along a line at the same angle as in the original geodesic data. So, for example, in illustrative image the red line is at on a bearing of of about 10 degrees relative to the red dot. That same angle is adopted in the projected image.

(5) To locate the blue dot, a ratio is determined between the sum of geodesic distances of the two green lines and the distance of the red line. Then, in the projected image, the blue dot is placed so that this ratio is the same. Note that the green lines will not be at 120-degrees from the red line in the projected image.

My question is whether this is a reasonable strategy for doing this kind of projection, and how I would go about implementing it in Mathematica.

• You can use this: mathematica.stackexchange.com/q/59463/5478
– Kuba
Jul 15, 2018 at 7:48
• Show us the code you have already tried, your definitions for maps and shapes and so on. This site is not about algorithms, but about Wolfram Mathematica programming so please edit your question to make it explicitly about programming in Mathematica/Wolfram Language . Include a formatted minimum example of the code you are working on. Jul 15, 2018 at 9:32
• @Kuba That is interesting, but it is not appropriate for a map transformation because the points move proportionately to the length of the perimeter which results in distortion. For example, in the case of the United States, the area around the great lakes would become gigantic in the projected map due to the longer coastline in that area. Jul 15, 2018 at 12:42
• @TylerDurden it is hard to avoid distortions here, isn't it? And you only mentioned your preferences somewhere in comments below answers. But yes, I agree it is only a quick start tip rather than full solution.
– Kuba
Jul 16, 2018 at 7:19

If the domain $\varOmega$ of the county is simply connect, one might use the Riemannian mapping theorem.

For $z_0 \in \varOmega^\circ$, we make the ansatz for the holomorphic map $f \colon \varOmega \to D^2$

$$f(z) = (z-z_0) \, \operatorname{e}^{u(z) + \operatorname{i}\!v(z)}.$$

Then $|f(z)| = |z-z_0| \, \operatorname{e}^{u(z)}$ has to equal $1$ for all $z \in \partial \varOmega$. The Cauchy-Riemann equations imply that $u$ is the unique solution to the following Poisson problem:

$$\begin{cases} \Delta u(z) &= 0, &z \in \varOmega^\circ,\\ u(z) &= - \log( |z-z_0|), &z \in \partial \varOmega. \end{cases}$$

# Implementation of $u$

This can be easily solved with the finite element method. Here is how this can be done (I chose a different country because the boundary curve of the USA seemed to have self-intersections):

First, we discretize the the interior of the domain:

Needs["NDSolveFEM"];
p = Most@Cases[CountryData["Germany", "Shape"], _Polygon, \[Infinity]][[1, 1, 1]];
R = ToElementMesh[
BoundaryMeshRegion[p, Line[Partition[Range[Length[p]], 2, 1, 1]]],
"MeshOrder" -> 1,
MaxCellMeasure -> {1 -> .01}
];
R["Wireframe"]


(*Initialization of Finite Element Method*)

vd = NDSolveVariableData[{"DependentVariables", "Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData[{"Space"} -> {R}];
cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{-IdentityMatrix[2]}},
"MassCoefficients" -> {{1}}
];
bcdata = InitializeBoundaryConditions[vd, sd, {DirichletCondition[u[x, y] == 0., True]}];
mdata = InitializePDEMethodData[vd, sd];

(*Discretization*)
dpde = DiscretizePDE[cdata, mdata, sd];
dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];
{load, stiffness, damping, mass} = dpde["All"];

(*Preparation of Dirichlet boundary conditions for u*)

bndedges = R["BoundaryElements"][[1, 1]];
bndvertices = Sort@DeleteDuplicates[Flatten[bndedges]];
bpts = R["Coordinates"][[bndvertices]];
z0 = {0., 0.};
b = ConstantArray[0., Length[R["Coordinates"]]];
b[[bndvertices]] = -0.5 Log[Total[(bpts - ConstantArray[z0, Length[bpts]])^2, {2}]];

(*Solving the system and creating an interpolating function*)

solver = LinearSolve[stiffness, Method -> "Pardiso"];
uvals = solver[b];
ufun = ElementMeshInterpolation[{R}, uvals];


# Implementation of $v$

Due to the Cauchy-Riemann equations, we have

$$\operatorname{grad}(v) = J \, \operatorname{grad}(u),$$

where $J$ is the rotation of 90 degrees in counter clockwise orientation. Since $\varOmega$ is simply connected, up to a constant, $v$ is defined uniquely by this relation. Hence, $v$ is also harmonic and is a solution of the following Neumann problem:

$$\begin{cases} \Delta v(z) &= 0, &z \in \varOmega^\circ,\\ \frac{\partial v}{\partial \nu} (z) &= \langle \nu(z) , J \operatorname{grad}(u)(z) \rangle, &z \in \partial \varOmega. \end{cases}$$ This pde is also amenable to a treatment with finite elements, although the preparation of the Neumann data needs a bit of fiddling. Notice that this second equation determines $v$ only up to a constant shift, which corresponds to the fact that $f$ is also unique up to a rotation of the disk.

(*Preparation of Neumann boundary conditions for v*)

gradu = {x, y} \[Function] Evaluate[D[ufun[x, y], {{x, y}, 1}]];
J = N@RotationMatrix[Pi/2];
p = R["Coordinates"];
{i, j} = Transpose[R["BoundaryElements"][[1, 1]]];
boundaryedgelengts = Sqrt[Total[(p[[i]] - p[[j]])^2, {2}]];
{α, β} = Transpose[bndedges];
vertexbndedgeconnectivity = SparseArray[
Transpose[{
Join[α, β],
Join[Range[Length[α]], Range[Length[β]]]
}] -> 1,
{Length[p], Length[bndedges]}
];

(*Solving the system and creating an interpolating function*)

b = vertexbndedgeconnectivity.(normalprojections boundaryedgelengts);
vvals = solver[b];
vfun = ElementMeshInterpolation[{R}, vvals];


Finally, we can obtain our approximation to the Riemannian mapping $f \colon \varOmega \to D^2$ by

f = {x, y} \[Function] Evaluate[
ComplexExpand[
ReIm[((x + I y) - (z0[[1]] + I z0[[2]])) Exp[ufun[x, y]] (Cos[vfun[x, y]] + I  Sin[vfun[x, y]])]
]
];


# Visualization

First, my favorite texture.

Now, the actual plots, using basically $f$ as texture map for $\varOmega$. (Plots were obtained with resolution MaxCellMeasure -> {1 -> .001}.)

tex = Import["https://i.stack.imgur.com/gRwc1.png"];
texcoords = Transpose[{
Total[(f @@@ p)^2, {2}],
ConstantArray[0.5, Length[p]]
}];
g = GraphicsRow[{
Graphics[{Texture[tex],
ElementMeshToGraphicsComplex[R, VertexTextureCoordinates -> texcoords]
}],
Graphics[{Texture[tex],
ElementMeshToGraphicsComplex[R, VertexTextureCoordinates -> texcoords,
"CoordinateConversion" -> (f @@@ # &)]
}]
},
ImageSize -> Full
]


• This is close to how i pelting uv maps was done in 3D in the past. Triangulate the mesh then make every edge a spring with rest length at starting length, then force the boundaries to the shape you want (circle) then let the spring dynamics handle the interiors Jul 15, 2018 at 14:31
• @joojaa I just wanted to brag around with my classical education ;) Indeed, these things are standard tasks in the domain of texture mapping -- which does not mean that they were easy. The art is to guarantee that there are no flipped triangles. There is much more to it than it would fit into a single post. Jul 15, 2018 at 14:34
• Nice use of the low level FEM functions. You might be interested in bndvertices === dbc["DirichletRows"] and normalprojections === MapThreadDot[ R["BoundaryNormals"][[1]], (gradu @@@ (0.5 (p[[i]] + p[[j]]))).(-J)]. The first part can be done in NDSolve directly, I'd need think about a way to compute the normal for the second part to be done in NDSolve. Jul 16, 2018 at 6:33
• You can also use something like: g = GraphicsRow[{Graphics[{Texture[tex], ElementMeshToGraphicsComplex[R, VertexTextureCoordinates -> texcoords]}], Graphics[{Texture[tex], ElementMeshToGraphicsComplex[R, VertexTextureCoordinates -> texcoords, "CoordinateConversion" -> (f[Sequence @@ ##] & /@ # &)]}]}, ImageSize -> Full] Jul 16, 2018 at 6:44
• @user21 Thank you so much! These tiny give-aways are actually really rewarding. The thing with the normal in the second equation was indeed bugging me for quite some time. Suggestions to simplify that are welcome. Jul 16, 2018 at 7:00

You could use a conformal map based on the Koebe–Andreev–Thurston circle packing theorem.

Image from Wikipedia.
It would not be a simple implementation from scratch, but it is very general, mapping any shape to any other. Used in computer graphics quite a bit, e.g., for texture mapping.

Image from Peter Schröder.
Maybe visit Ken Stephenson's circle packing web pages.

• How is this implemented in Mathematica? Jul 15, 2018 at 14:06
• @TylerDurden: I only know it has been implemented (perhaps not in Mathematica). I added a link. Jul 15, 2018 at 14:15
• The Stephenson suffers from the same problem mentioned in my comment to the original question in that it distorts the map in areas of complex perimeters. The "great lakes problem". Eg in the example shown you can see how the gray dot is distorted. In a map projection this point should be relatively close to the coastline, because it is close to the coastline in the geodesic data. You may want to read the question more closely and think about how maps are designed to answer better.. Jul 15, 2018 at 14:35
• @Tyler Any such mapping would introduce distortion of distances. Conformal maps have the nice property that they still preserve angles. Jul 15, 2018 at 14:40
• @HenrikSchumacher For my purposes, I care most that the distance to the nearest point on the coastline are close to their geodesic distance. The angle to a distant point is not important. Jul 15, 2018 at 14:43