Draw a multi polygon with metric dimensions over part of a sphere, using d3-geo

Question

I have the following need:

given a multi-polygon with metric dimensions, I need to draw it over a specific part of a sphere and then output it in GeoJSON format.

additionally, the "part of the sphere" is a rectangle made of 4 points, which could be rotated in relation to the Equator line.

as an example, imagine a MultiPolygon with a rectangle in it, like this, where the numbers are in meters:

[
  [
    [
      [0, 10],
      [0, 0],
      [20, 0],
      [20, 10]
    ]
  ]
]

as you can see, such polygon is aligned in right angle with the axis.

Then, given any rectangular area of the sphere, for instance the one below:

{
    leftTop: [ 13.377948, 52.521293 ],
    leftBottom: [ 13.377962, 52.521250 ],
    rightBottom: [ 13.378174, 52.521275 ],
    rightTop: [ 13.378161, 52.521318 ]
}

I want to translate, rotate and fit that multi-polygon in this area.

So, d3-geo package seems to offer the proper tools, such as projection with the combination of scale/translate/rotate/fitExtent with invert and so on. But I can't find the right combination to deal with rotated rectangular areas.

As the vast majority of examples of d3-geo in the web are either from GeoJSON to pixels or using d3 as an SVG rendering engine, etc. it's being hard for me to abstract my needs into function calls.

Does anyone have a light to help me solving this? I guess there must be a proper term in Geometry/Trigonometry to describe such calculation.

I have used for my self "projecting a plane polygon over a geodesic area", but I guess there must be a better and more specific term. At least if I find the such term it would be a good start.

Update:

I made up some illustrations of what I'm trying to achieve. Keep in mind that even though they are visual images here, in my case I'm not looking for visual rendering tools, as I will do just a transformation from Array/Objects into GeoJSON.

1. A polygon I want to project. Notice that the light gray box is showing the virtual area where such polygon is positioned over, so, as the polygon is in the middle, it should be projected into the middle of the geographic area too.

2. A rectangular sector over a geoid (the Earth in this case). Notice the red labels pointing the reference corners, which should indicate the rotation.

3. The result, with the polygon projected into the given area (which would be in fact in GeoJSON format, with respective coordinates).

Answer 1

If I understand the problem correctly, I believe I have a solution. However, the question is quite broad without any code, consequently, my answer will be fairly broad, but should clearly outline an approach to achieve the desired result.

Based on the images the problem appears to be: How to fit shapes into geographic bounding boxes where the bounding box is projected into planar space and may have a rotation.

If this is accurate, then my answer should hopefully be useful. For sake of clarity, when I use the word shape I'm refering to the non-geographic polygon.

The general pattern of the solution is:

1. Get the orientation of the reference polygon
2. Rotate the map - not the shape
3. Using a geoIdentity, fit the shape to the geographic bounding box

4. a. Finish by applying the opposite rotation as you started with,

   b. Rotate the entire svg/canvas to counter the initial rotation; or,

   c. Keep the rotated map.

---

Getting the Orientation of the Geographic Rectangle

You need to establish the rotation of the geographic feature relative to the projection plane, this appears to be the crux of the issue (You state the rotated rectangular area of a sphere, but areas can only be rectangular in projected space, and based off the images, this interpretation appears to be correct).

D3 can be a bit challenging in this regard as it does not render paths along Cartesian coordinates, but interpolates great circle distances. This will be more of an issue with country or continent sized geographic features, but should be negligible at city scale.

If starting with the top left point ([x1,y1]) and using the bottom right point ([x2,x1]), you should be able to determine the angle of this line in relation to what it should be: a vertical line starting at the first coordinate and ending at the second. Given that winding order matters with d3, if you have one point, the second coordinate will always correspond to the same vertex of the rectangle.

The method to get this angle is fairly straightforward:

var p1 = [long,lat];   // geographic coordinate space for the two points
var p2 = [long,lat];

var x1 = projection(p1)[0];  // projected svg coordinate space for the two points
var x2 = projection(p2)[0];
var y1 = projection(p1)[1];
var y2 = projection(p2)[1];

var dx = x1 - x2;  // run
var dy = y1 - y2;  // rise

var angle = Math.atan(dx/dy);  // in radians, multiply by 180/π to get degrees

We're projecting the two coordinates, calculating the projected differences in x and y coordinates, measured in pixels. Run it through Math.atan() and we have an angle.

But wait, there's more.

Rotating the Map

The method we've used to calculate the angle is fine, but we need to modify it to rotate the map. If we rotate the map, it will rotate around [0,0] (long,lat), the default rotational center of most projections. We need to center the map at [-x,-y] where x and y represent the mid point of the two points used to calculate orientation, measured in geographic coordinate space (not projected).

We need to center the map by rotation before we calculate angle, as this will alter the angle. For this we need d3.geoInterpolate which calculates points from p1 (0) to p2 (1), we need the half way point, so we feed it 0.5:

var pMid = d3.geoInterpolate(p1, p2)(0.5);

Now we can apply it to the projection prior to calculating angle:

projection.rotate([-pMid[0],-pMid[1]]);

Why negative values? We move the earth under us

Now we can go through the angle calculation. Once we have the angle, we can apply it:

projection.rotate([-pMid[0],-pMid[1],-angle])  // angle in degrees

We're half way there, using the image below, we used the geographic coordinates of A and B to determine the geographic center C. Then using the projected coordinates of A and B we determine angle α which we then use to to rotate the projected coordinates so that line AB is vertical on the map.

Fit the Shape to the Geographic Bounding Box

So we have half the problem solved, now we need to project the shape. We will use a second projection for the shape, a plain geoIdentity would work, it allows us to use a fitSize or fitExtent method while projecting coordinates with no transform. Just note you probably want to flip this feature on the y axis: svg y values start with 0 at the top, more standard Cartesian y values start at the bottom.

We'll want to use fitExtent, which will allow us to set a bounding rectangle for the shape. proejction.fitExtent([[x,y],[x,y]],feature) takes an array containing the upper left and the bottom right corners of a bounding box (in svg coordinates) to hold a feature (geojson feature).

Keep in mind that projected geographic feature we righted will be more rectangular the smaller it is in geographic (not projected) size, larger areas might have less square properties, especially in relation to the right hand side of the rectangle. For large geographic areas you might need to revise the rotation calculation but at some point projected rectangles just don't align with "rectangles" on a sphere.

To get the bounding box for fitExtent we can use path.bounds():

Returns the projected planar bounding box (typically in pixels) for the specified GeoJSON object. The bounding box is represented by a two-dimensional array: [[x₀, y₀], [x₁, y₁]], where x₀ is the minimum x-coordinate, y₀ is the minimum y-coordinate, x₁ is maximum x-coordinate, and y₁ is the maximum y-coordinate. (API docs)

Great, now we have two bounding boxes that should have the same points, we use:

projection2.fitExtent(path.bounds(geoRectangle));

Now we have overlain the shape on the geographic feature:

var feature = { "type": "Feature","properties": {},"geometry": {"type": "Polygon","coordinates": [ [ [-81.62841796875,24.307053283225915],[-75.9375,21.88188980762927],[ -77.3876953125,18.8543103618898],[       -83.1884765625,21.268899719967695 ], [-81.62841796875,24.307053283225915]]]}};
var triangle = {"type": "Polygon","coordinates": [[[0, 0], [10, 10], [20, 0], [0, 0]]]};


var width = 500; var height = 300;
var svg = d3.select("body")
  .append("svg")
  .attr("width",width)
  .attr("height",height);
  
var projection = d3.geoMercator().scale(500/Math.PI/2).translate([width/2,height/2]);
var path = d3.geoPath().projection(projection);


d3.json("https://unpkg.com/world-atlas@1/world/110m.json", function(error, world) {
  if (error) throw error;

  var p1 = feature.geometry.coordinates[0][1]; // first point in geojson
  var p2 = feature.geometry.coordinates[0][2]; // second point in geojson
  
  var pMid = d3.geoInterpolate(p1, p2)(0.5);   // halfway point between both points
  
  projection.rotate([-pMid[0],-pMid[1]]);                          // rotate the projection to center on the mid point
  projection.fitExtent([[135,135],[width-135,height-135]],feature) // optional: scale the projected feature, may offer benefits for very small features
 
  var dx = projection(p1)[0] - projection(p2)[0]; // run, difference between projected points x values
  var dy = projection(p1)[1] - projection(p2)[1]; // rise, difference between projected points y values
  
  var a = Math.atan(dx/dy) * 180 / Math.PI; // get angle and convert to degrees
  
   projection.rotate([-pMid[0],-pMid[1],-a]);  // adjust rotation to straighten feature
   projection.fitExtent([[135,135],[width-135,height-135]],feature) // scale and translate the feature.
   
  // draw world map, draw feature
  svg.append("path")
    .attr("d",path(topojson.mesh(world)))
    .attr("fill","none")
    .attr("stroke","black")
  
  svg.append("path")
    .attr("d", path(feature))
    .attr("fill","none")
    .attr("stroke","steelblue");


  // set up the projection and path for the shape       
  var projection2 = d3.geoIdentity().reflectY(true);
  var path2 = d3.geoPath().projection(projection2);

  // scale the shape's projection for the shape using the bounds of the geographic feature
  projection2.fitExtent(path.bounds(feature),triangle);

  // draw the shape
  svg.append("path")
    .attr("d", path2(triangle));
  
});

<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/4.11.0/d3.min.js"></script>
<script src="https://unpkg.com/topojson-client@3"></script>

This is a hand drawn polygon, I placed a business card on the monitor and traced it, hence the greater than normal imperfections

End Game

Now you have a shape drawn and overlain on top of some geographic feature, hopefully fairly decently. Now what? Technically this solves the key issue, but we've created a new one if you don't want a tilted map. Options for solving this are rotating the entire svg/canvas or taking each point in the shape and reproject them and turning them into geographic coordinates rather than boring old Cartesian ones. There are others, but these two seem most straightforward.

If dealing with images, you could reproject the image pixel by pixel, but don't expect this to be quick, see this block. This answer looks at images and fitting them to a projection if they have known bounds, with very small distances a Mercator should work fine as the projection.

Reprojecting points and GeoIdentity

As indicated in the comments, d3.geoIdentity makes reprojection impossible if following the pattern: var latlong = projection.invert(projection2([x,y])), as a geoIdentity doesn't return a function, just an object. Instead, we can define the behavhior (flipping on the y, and using identity scale and transform):

projection2.project = function([x,y]) {
  var s = this.scale();
  var t = this.translate();
  return [(x * s) + t[0] , -y * s + t[1]]
}

Which in cursory testing appears to return the projected point, which is the standard behavior of a geoProjection not included in a geoIdentity. The usage pattern should look like:

projection.invert(projection2.project([x,y]))

Answer 2

I found a way to implement what I need, although I didn't use d3-geo, but used instead a combination of Turf.js and geodesy libraries.

Here is it (search for "projectPlainOverGeoArea"):

https://gist.github.com/marinho/8a7e71b53c7da97f9579146742bb54de

The code does the following:

1. calculate width and height in meter distances from the geographical rectangle (the area I want to project to). It uses geodesy spherical/haversine functions

2. project each point of each path of the polygon onto the equivalent geo coordinates by using the proportional distances from the polygon rectangle onto the geo rectangle. The result here is the polygon the expected measures aligned with the Equator line.

3. uses Turf's transformRotate to rotate the projected polygon using left/bottom as pivot

That seems to work fine in the tests I've made, but to be honest I'm unhappy about it, as I feel there must a more elegant solution for it with d3-geo or another library.