D3 Geojson Geocircle Ellipse Equivalent

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 GeoJS

Solution 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:

functioncreateEllipse(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 ellipsesvar coords = [];

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

    // get the radius at that anglevar 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);


functioncreateEllipse(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 anglesvar coords = [];
	
	for (var i = 0; i <= n; i++) {
	    // get the current anglevar θ = Math.PI*2/n*i + rotation;
		
		// get the radius at that anglevar 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] };
	
}

<scriptsrc="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:

functioncreateEllipse(a,b,x=0,y=0,rotation=0) {
	
	var k = Math.ceil(36 * (Math.max(a/b,b/a))); // sample anglesvar coords = [];
	
	for (var i = 0; i <= k; i++) {
	
		// get the current anglevar angle = Math.PI*2 / k * i + rotation
		
		// get the radius at that anglevar 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] };
}
 
functiongetLatLong(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) { returnprojection(d)[0] })
  .attr("cy", function(d) { returnprojection(d)[1] })
  .attr("r", 4)
  .attr("fill","orange");

<scriptsrc="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