d3 GeoJSON geoCircle ellipse equivalent

Question

The title pretty much says it all. I'm looking for a convenient way to generate a geoJSON polygon defining an ellipse similar to d3-geo's d3.geoCircle()(); I want to use this GeoJSON ellipse with d3-geo. To clarify with and example, Cesium has this capability with a simple function allowing you to create an ellipse like so:

var ellipse = new Cesium.EllipseGeometry({
  center : Cesium.Cartesian3.fromDegrees(-75.59777, 40.03883),
  semiMajorAxis : 500000.0,
  semiMinorAxis : 300000.0,
  rotation : Cesium.Math.toRadians(60.0)
});

If that function returned GeoJSON I'd be set. What's the best way to generate a GeoJSON polygon defining an ellipse?

Answer 1

D3 doesn't offer anything that can really help here. Vanilla javascript can achieve this fairly easily. First let's create a geojson ellipse in Cartesian coordinate space. After, we can use the haversine formula to draw the ellipse.

1. Create a geojson ellipse in Cartesian coordinate space.

This is pretty straightforward, the method I'm using is to calculate the radius of the ellipse at a given angle. Using these polar coordinates we can stitch together an ellipse. The formula for the radius of an ellipse at a given point can be found pretty easily, I used this source, which gives us:

So, we can easily iterate through a series of angles, calculate the radius at that angle, and then translate this polar coordinate into a Cartesian coordinate. Perhaps something like:

function createEllipse(a,b,x=0,y=0,rotation=0) {

  rotation = rotation / 180 * Math.PI;
  var n = n = Math.ceil(36 * (Math.max(a/b,b/a))); // n sampling angles, more for more elongated ellipses
  var coords = [];

  for (var i = 0; i <= n; i++) {
    // get the current angle
    var θ = Math.PI*2/n*i + rotation;

    // get the radius at that angle
    var r = a * b / Math.sqrt(a*a*Math.sin(θ)*Math.sin(θ) + b*b*Math.cos(θ)*Math.cos(θ));

    // get the x,y coordinate that marks the ellipse at this angle
    x1 = x + Math.cos(θ-rotation) * r;
    y1 = y + Math.sin(θ-rotation) * r;

    coords.push([x1,y1]);
  }

  // return a geojson object:
  return { "type":"Polygon", "coordinates":[coords] };

}

Note: a/b: axes (in pixels), x/y: center (in pixels), rotation: rotation in degrees

Here's that in a quick snippet:

var geojson = createEllipse(250,50,200,200,45);

var svg = d3.select("body")
  .append("svg")
  .attr("width",600)
  .attr("height",500);
  
var path = d3.geoPath();

svg.append("path")
 .datum(geojson)
 .attr("d",path);


function createEllipse(a,b,x=0,y=0,rotation=0) {

 rotation = rotation / 180 * Math.PI;
 var n = n = Math.ceil(36 * (Math.max(a/b,b/a))); // n sample angles
 var coords = [];
 
 for (var i = 0; i <= n; i++) {
     // get the current angle
  var θ = Math.PI*2/n*i + rotation;
  
  // get the radius at that angle
  var r = a * b / Math.sqrt(a*a*Math.sin(θ)*Math.sin(θ) + b*b*Math.cos(θ)*Math.cos(θ));
  
  // get the x,y coordinate that marks the ellipse at this angle
  x1 = x + Math.cos(θ-rotation) * r;
  y1 = y + Math.sin(θ-rotation) * r;

  coords.push([x1,y1]);
 }
 
 // return a geojson object:
 return { "type":"Polygon", "coordinates":[coords] };
 
}

<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/4.13.0/d3.min.js"></script>

2. Apply the haversine formula.

One of the best resources on the haversine and related functions I know of is at Moveable Type Scripts. The formula I have came from there a few years back and has had a few cosmetic modifications. I'm not going to break down the formula here, as the linked reference should be useful.

So, rather than calculating the Cartesian coordinates, we can take the polar coordinate and use the angle as bearing and the radius as distance in the haversine formula, which should be relatively trivial.

This could look like:

function createEllipse(a,b,x=0,y=0,rotation=0) {
 
 var k = Math.ceil(36 * (Math.max(a/b,b/a))); // sample angles
 var coords = [];
 
 for (var i = 0; i <= k; i++) {
 
  // get the current angle
  var angle = Math.PI*2 / k * i + rotation
  
  // get the radius at that angle
  var r = a * b / Math.sqrt(a*a*Math.sin(angle)*Math.sin(angle) + b*b*Math.cos(angle)*Math.cos(angle));

  coords.push(getLatLong([x,y],angle,r));
 }
 return { "type":"Polygon", "coordinates":[coords] };
}
 
function getLatLong(center,angle,radius) {
 
 var rEarth = 6371000; // meters
 
 x0 = center[0] * Math.PI / 180; // convert to radians.
 y0 = center[1] * Math.PI / 180;
 
 var y1 = Math.asin( Math.sin(y0)*Math.cos(radius/rEarth) + Math.cos(y0)*Math.sin(radius/rEarth)*Math.cos(angle) );
 var x1 = x0 + Math.atan2(Math.sin(angle)*Math.sin(radius/rEarth)*Math.cos(y0), Math.cos(radius/rEarth)-Math.sin(y0)*Math.sin(y1));
 
 y1 = y1 * 180 / Math.PI;
 x1 = x1 * 180 / Math.PI;
   
 return [x1,y1];
} 

// Create & Render the geojson:
var geojson = createEllipse(500000,1000000,50,70); // a,b in meters, x,y, rotation in degrees.
var geojson2 = createEllipse(500000,1000000)

var svg = d3.select("body")
  .append("svg")
  .attr("width",600)
  .attr("height",400);
  
var g = svg.append("g");

var projection = d3.geoMercator().translate([300,200]).scale(600/Math.PI/2);

var path = d3.geoPath().projection(projection);

g.selectAll("path")
 .data([geojson,geojson2])
 .enter().append("path")
 .attr("d", path);
 
g.selectAll("circle")
  .data([[50,70],[0,0]])
  .enter().append("circle")
  .attr("cx", function(d) { return projection(d)[0] })
  .attr("cy", function(d) { return projection(d)[1] })
  .attr("r", 4)
  .attr("fill","orange");

<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/4.13.0/d3.min.js"></script>

Note: a/b axes in meters, x,y,rotation in degrees

That's a pretty boring demonstration, perhaps this simple demonstration is better:

The formula I'm using assumes a earth is a sphere, not an ellipsoid, this can lead to errors in distance of up to 0.3%. However, depending on map scale, this will often be less than the stroke width.

I might have to try and make a particularly visually challenging version of a tissot's indicatrix with this

Snippets use default parameter values that are not compatible with IE, example block offers IE support