Nord de Guerre
Here’s a survey plat made in 1919 by a US Army unit in France.
In the aftermath of WW1, teams of American military surveyors produced these as they went around France mapping the grave sites of fallen soldiers.
On this plat we can see a number of bearings and distances from the I.P. (Initial Point), which is a small triangle at the centre of the image above. The bearing to the grave of Private Joseph A. Cain, for example, is “S 55°40′ W 598.5.” This means 55° 40′ west of due south — although we have to bear in mind that it was magnetic south, and specifically the magnetic south of 1919 in this location — at a distance of 598.5 feet.
The “I.P.” itself is a copper nail whose approximate coordinates tie the whole plat to a larger coordinate system. But its coordinates are given as the enigmatic N 296.80, E 309.02. What do these numbers mean?
They are French Lambert coordinates of the World War 1 era, more properly known as coordinates in the Nord de Guerre projection. They were used extensively in the mapping done by the French, British and American armies during the war. Clifford Mugnier, in his 2001 article on the history of grids and datums in the French Republic, gives us the parameters of the projection:
On the other hand, the French Nord de Guerre Zone (1914-1948) was used, and it was based on the French Army Truncated Cubic Conic where the Latitude of Origin was φ0 = 49° 30′ 00″, the Central Meridian was λ0 = 7° 44′ 13.95″ East of Greenwich, the Scale Factor at Origin (m0) = 0.999509082, and the False Easting was 500 km and the and False Northing was 300 km. The ellipsoid of reference was the Plessis Reconstituted.
Mugnier, Clifford, The French Republic, PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING, January 2001, pp. 33-35
These are welcome numbers to anyone who might want to re-create this projection. But if you aren’t doing that, you could still summarize it as follows: Lambert coordinates are in kilometers, and they are relative to an origin point at 49° 30′ N, 7° 44′ 13.95″ E. However, this origin is not N 0.00,E 0.00, but has been assigned N 300.00, E 500.00. This ensures that all coordinates within the area they are mapping will be positive.
(Ironically, the origin of this coordinate system was inside Germany.)
If having a central meridian at 7° 44′ 13.95″ east of Greenwich seems weird to you, understand that the French didn’t use the Greenwich Meridian. They used the Paris Meridian — a line through the Paris Observatory — which was 2° 20′ 13.95″ (or 2.33720833333333°) east of Greenwich. The central meridian of this Lambert coordinate system was then a further 5° 24′ east of that. That’s not a logical number either, until you realize that the French also weren’t using degrees: they were using grads. A grad is one four-hundredth of a circle, so 5° 24′ is exactly 6 grads. Hence the central meridian is the common-sense “six grads west of Paris.” If you add 5° 24’ to 2° 20′ 13.95″ you get 7° 44′ 13.95”.
In QGIS, you can see that the PROJ library comes with a definition for this historic projection: “ATF (Paris) / Nord de Guerre.” It looks like this:
+proj=lcc +lat_1=49.5 +lat_0=49.5 +lon_0=5.4 +k_0=0.99950908 +x_0=500000 +y_0=300000 +a=6376523 +rf=308.64 +pm=2.33720833333333 +units=m +no_defs
In that PROJ definition you can see the parameters that Mugnier described. lat_1 and lat_0 are 49.5°, a.k.a. 49° 30′ . Lon_0 is only 5.4° but it’s added to a prime meridian (+pm) of 2.33720833333333 (Paris Meridian), to make the 7° 44′ 13.95″ E. The scale factor (+k) is 0.99950908, and the false eastings and northings (+x_0 and +y_0) are given. The +a=6376523 and +rf=308.64 are the parameters of the Plessis Reconstituted ellipsoid.
But this is not the problem.
Can we transform?
Last year, I was contacted by a database administrator volunteering for a non-profit in the United States called DoughboyMIA. They attempt to locate the remains of U.S. soldiers who died in France in WW1. In many cases these remains were disinterred in the 1920s and transferred to a French military cemetery, or even to a cemetery in the United States. (Private Joseph Cain is one of these: his remains were moved to a cemetery in Leadville, Colorado.) However, many others remain buried where the plats show them.
Would it be possible, my friend asked, to convert these old Lambert coordinates into lat/long, so they could call them up in Google Earth and then figure out where the burial sites were? They had been trying to do this using an overlay grid on smaller scale battlefield maps, but they were generally getting no closer than 100 metres.
I thought it sounded pretty simple. It was not.
I thought I would just use QGIS, and feed in the Nord de Guerre coordinates on one end and get out WGS84 lat/longs at the other. Modern GIS software does this all the time: we assemble data layers in various Coordinate Reference Systems, and it automatically converts them so everything is in the CRS of your map project.
The PROJ package, which underlies QGIS’s ability to convert coordinate systems, contains a command-line utility called cs2cs (“coordinates system to coordinate system”) that does this sort of thing. To run it, we will need the EPSG codes for Nord de Guerre and WGS84, which are 27500 and 4326, respectively. (You can look these up in the list of coordinate reference systems in QGIS.) Then, we have to feed in the Lambert coordinates in metres, so, grabbing N 296.80, E 309.02 from the plat above, we move the decimal point three places to the right: 296.80 becomes 296800 and 309.02 becomes 309020. We will present these to cs2cs in the conventional order: easting (x-axis) first, northing (y-axis) second.
$ echo 309020 296800 | cs2cs -f "%.6f" EPSG:27500 +to EPSG:4326
49.441146 5.101183 0.000000
We get out latitude 49.441146° North, longitude 5.101183° East. (The switch -f “%.6f meant format the answer to six decimal places.)
That was easy. However, if you paste “49.441146 5.101183” into the search bar on Google Maps, you get both good news and bad news.
The good news is that the place is recognizably the area surveyed in Plat A-6. You can see the farmhouse, the road, the intersection. After 100 years, it’s all still there. We’re definitely in the ballpark. The bad news is that the location of the I.P. is off a lot. It’s supposed to be right at a junction of two roads. But cs2cs put it 140 metres to the northeast.
As I was hunting around for why, I discovered that one can ask PROJ about how it does a transformation between any pair of coordinate systems. In this case (27500 to 4326), it gave this unsettling answer:
$ projinfo -s EPSG:27500 -t EPSG:4326
Candidate operations found: 1
Operation No. 1:
unknown id, Inverse of Nord de Guerre + Transformation from ATF (Paris) to ATF (Paris) altered to use prime meridian of WGS 84 + Ballpark geographic offset from ATF (Paris) altered to use prime meridian of WGS 84 to WGS 84, unknown accuracy, World, has ballpark transformation
The bad news here is “ballpark.” Ballpark means no one has defined how to transform from Nord de Guerre coordinates to WGS84 lat/long. Specifically, the problem is moving from the ellipsoid and datum for Nord de Guerre to the ellipsoid and datum for WGS84. PROJ is forced to assume that Nord de Guerre was based on the WGS84 ellipsoid, which it wasn’t. So the result is 140 metres out.
To make sense of this, I just have to digress a bit into the history of datums and ellipsoids. Right up through the mid 20th century, mapping all over the world used whatever ellipsoids and datums the local national mapping program preferred. There were many ellipsoids (e.g., Bessel, Krassowsky, Clarke-1866) and many datums (e.g., Djakarta, Pulkovo-1942, Old Hawaiian). Each ellipsoid was an almost-spherical model of the earth, defined by its major axis and its minor axis. Each datum was, in effect, the location of 0° N, 0° E on that ellipsoid.
But with the advent of GPS, the idea arose to create some kind of overarching world ellipsoid and datum, and then figure out how to transform the coordinates of the many local datums and ellipsoids into that common-ground datum. That world ellipsoid/datum was WGS84, and you can see all the transformations that were figured out to link it to the pre-existing datums and ellipsoids in, among other places, a document from the International Hydrographic Association.
For example, if you wanted to transform coordinates in the Old Hawaiian datum, which was based on the Clarke 1866 ellipsoid, it was all worked out how you get them to WGS84 on page 79.
The transformation parameters of the “mean solution” are ΔX=61, ΔY=-285 and ΔZ=-181. (I’ll explain what that means below.)
This knowledge is baked into PROJ. If you look up the Old Hawaiian coordinate system in QGIS you’ll find its EPSG number is 4135, and the PROJ string looks like this:
+proj=longlat +ellps=clrk66 +towgs84=61,-285,-181,0,0,0,0 +no_defs
See the +towgs84 part? There are the ΔX, ΔY and ΔZ from the chart above: 61,-285 and -181. They tell PROJ how to transform coordinates from Old Hawaiian CRS into WGS84. So if you happen to have some data projected in the Old Hawaiian projection, you can display it over a GPS track recorded yesterday (i.e., WGS84 coordinates) and everything should line up.
Sadly, though, no one appears to have figured out the ΔX, ΔY and ΔZ for coordinates coming from the Plessis Reconstituted ellipsoid, and specifically not one for the Guerre de Nord projection. Both the ellipsoid and the projection had just been out of use too long by the time WGS84 came around. (France scrapped Guerre de Nord right after WW1 and came up with a completely new Lambert system, based on the Clarke 1880 ellipsoid.)
How do these transformations even work?
In the +towgs84 expression above, you can see that there is a list of seven parameters, but only the first three are used (“61,-285,-181,0,0,0,0”). These are the parameters of a Helmert Transformation, and most Helmert transformations that were calculated for getting old coordinate systems to WGS84 use only the first three parameters, called ΔX, ΔY and ΔZ. You can call these transformations “3-parameter Helmerts.”
I’ll explain how they are used (which I learned largely from the PROJ page about geodetic transformations), but first you have to picture the ECEF (Earth Centered Earth Fixed) coordinate space.
The ECEF is a kind of common ground where ellipsoids meet and differences can be quantified. ECEF, which no one would ever use for mapping or navigation, is a Cartesian coordinate system, like the graphs with X- and Y-axes you made in high school. Being in three dimensions, it has three mutually perpendicular axes, X, Y and Z, that meet at the center of the Earth.
The Z-axis runs straight up through the North Pole. The X-axis comes out through the point at 0° N, 0° E (in the Atlantic, just south of Accra, Ghana). The Y-axis comes out through the point at 0° N, 90° E (in the Indian Ocean between Sri Lanka and Sumatra).
ECEF is quite weird. We’re used to thinking of places being located in spherical coordinates, where the latitude is the angle up from the equator (φ in the diagram above) and the longitude (λ) is the angle east from the Greenwich meridian. For example, the Greenwich Observatory is at 51° 28′ 40″ N, 0° E.
But in ECEF coordinates, every place is located with its (X, Y, Z) coordinates, relative to the earth’s centre. The Greenwich Observatory is at (3,979,301, 0, 4,966,347) — which is like, “3,979,301 metres out on the x-axis, no metres on the y-axis, and 4,966,347 metres up on the z-axis.”
Why would anyone want this? Well, it’s a good place mathematically to step from one ellipsoid to another. And here’s how the 3-parameter Helmert is usually done, going from some old ellipsoid and datum to WGS84:
- You get a place’s latitude and longitude on its old ellipsoid.
- You convert that latitude and longitude into an (X, Y, Z) coordinate in ECEF. (You can do this because you have the equations that define the ellipsoid.)
- You apply the three Helmert parameters (ΔX, ΔY and ΔZ) by adding them to their respective coordinates. In the case of Old Hawaiian above, you would add 61 metres to the X coordinate, subtract 285 metres from the Y coordinate, and subtract 181 metres from the Z coordinate.
- Now you convert this (X, Y, Z) to a latitude/longitude on the WGS84 ellipsoid. (Again, you can do this because you have the equations that define the ellipsoid.)
You can see that if we want to go from Nord de Guerre to WGS84, we really have all of the information we need except the three Helmert parameters.
Let’s figure out the parameters
On the wikipedia page about Helmert Transformations I read that you could estimate your own parameters if you had a series of points whose coordinates in both reference systems were known. I realized this would not be too hard to arrange, because once I had the ballpark of where that plat A-6 was, it had been easy to georeference it to existing features on satellite imagery. I was able to see what the WGS84 coordinates of the I.P. at N 296.80, E 309.02 should be.
And the US National Archives has many of these survey plats from 1919 and 1920, from all over northeast France.
So let’s walk through converting the first pair of Lambert coordinates to the ECEF coordinates, and converting the corresponding WGS84 coordinates into ECEF as well. We can compare them and see what the differences are (X, Y and Z). Those differences will be our Helmert parameters, and if we repeat this process on many points we should be able to average our Helmert parameters to get a “mean solution.”
Step one is to convert N 296.80, E 309.02 to lat/long, and by this we mean lat/long on the Plessis ellipsoid. You can go to epsg.io and look up any projected coordinate reference system (CRS) and find the unprojected (geodetic) CRS that underlies it.
In this case it’s “ATF (Paris)” whose EPSG number is 4901. (ATF is the Ancienne Triangulation Francaise).
We can convert the Nord de Guerre easting and northing to ATF (Paris) lat/long using cs2cs like this:
$ echo 309020 296800 | cs2cs -f "%.6f" EPSG:27500 +to EPSG:4901
49.441146 2.763974 0.000000
(Notice that it gives us a third coordinate in the result, which is elevation above the ellipsoid. It’s zero and we’re going to ignore it.)
So now we’re on the Plessis ellipsoid. The longitude in the result, 2.763974° seems wrong, and that’s because we still have to add the Paris meridian (2.33720833333333) to it. And that gives us 5.101182°.
Step two is to convert 49.441146° N, 5.101182° E on the Plessis ellipsoid to XYZ in the ECEF, like this:
$ echo 49.441146 5.101182 | cs2cs EPSG:4901 +to +proj=cart +ellps=plessis
4137509.16 369348.89 4822160.05
That’s a funky set of numbers, eh? (4137509.16, 369348.89, 4822160.05)
Now let’s work backwards from the WGS84 coordinates that I pulled off the georeferenced plat, the coordinates we want to get: 49.44049352° N, 5.09992774° E. We convert to ECEF (step four, backwards) with
$ echo 49.440493520 5.099927740 | cs2cs EPSG:4326 +to '+proj=cart +ellps=WGS84'
4138889.00 369380.74 4822556.01
That was easy. Now we compare the coordinates. A spreadsheet is useful.
For all three coordinates, the values on the WGS84 ellipsoid are bigger, so our (ΔX, ΔY, ΔZ) are all positive: (1379.84, 31.85, 395.96).
To make a long story short, I repeated this with twelve more IP points on different plats scattered around northeast France.
In the process I learned a few things. Some plats, usually those showing actual walled graveyards, did not have an IP, but instead had coordinates that corresponded to something else, like the graveyard entrance or the church steeple. And most plats had a Lambert coordinate that was plainly rounded to the nearest hundred metres (like “N 296.80“). In other words, these guys didn’t have GPS — they were triangulating from benchmarks — and mostly they were only confident to within 100 metres. So on the whole I did not expect accuracy better than 100 metres.
To my surprise (and delight), the differences in ECEF were fairly tightly clustered.
ΔX was almost always in the high 1300s. ΔY was typically near 40, and had the highest standard deviation. ΔZ was generally quite close to 390. I could expect most of the time to be well within 100 metres.
When I averaged these, I got ΔX=1383.8, ΔY=38.7, and ΔZ=392.0.
My +towgs84 term for converting Nord de Guerre to WGS84 would then be:
+towgs84=1383.8,38.7,392,0,0,0,0
Employing the transformation
If you wanted a full PROJ definition for Nord de Guerre that you could add to QGIS as a custom projection, a projection which would know how to transform itself into WGS84, it would be
+proj=lcc +lat_1=49.5 +lat_0=49.5 +lon_0=5.4 +k_0=0.99950908 +x_0=500000 +y_0=300000 +a=6376523 +rf=308.64 +pm=2.33720833333333 +towgs84=1383.8,38.7,392,0,0,0,0 +units=m +no_defs
But we didn’t actually use that. DoughboyMIA doesn’t use QGIS. They just have long lists of coordinates. And for purposes of converting large batches of coordinates, it was more useful to provide a command line, and direct its input to come from a file .
For example you could have a file called “lamberts.csv” that looks like this. It’s pairs of Lambert coordinates separated by tabs: northing, easting, elevation. (I always set elevation to 0.)
309020 296800 0
300960 276980 0
452380 96610 0
327630 280480 0
361560 239430 0
299100 290000 0
297700 294410 0
307680 280000 0
327660 280510 0
294760 269960 0
237810 222600 0
404970 165890 0
446330 347930 0
At this point I had begun to prefer the command-line utility cct to cs2cs. cct is another tool from PROJ; the chief difference between it and cs2cs is that cct expects an elevation parameter as a third coordinate. The “< lamberts.csv” at the end is me passing it the tab-separated file of input Lambert coordinates.
cct -I +proj=lcc +lat_1=49.5 +lat_0=49.5 +lon_0=5.4 +k_0=0.99950908 +x_0=500000 +y_0=300000 +a=6376523 +rf=308.64 +pm=2.33720833333333 +towgs84=1383.8,38.7,392,0,0,0,0 +units=m +no_defs < lamberts.csv
And the WGS84 output looks like this. Longitude first, then latitude; then two coordinates which I discard: elevation and time.
4.9988362488 49.2596692455 0.2851 0.0000
7.1011763673 47.6678928315 -0.8421 0.0000
5.3637106750 49.2992738854 -0.1256 0.0000
5.8444804118 48.9387726823 -0.2154 0.0000
4.9667099639 49.3761059620 0.1611 0.0000
4.9451843243 49.4152856536 0.1278 0.0000
5.0896851527 49.2889797249 0.1609 0.0000
5.3641101466 49.2995521280 -0.1264 0.0000
4.9172892218 49.1945037089 0.4440 0.0000
4.1672182134 48.7471097765 1.5824 0.0000
6.4542969576 48.2855823021 -0.3381 0.0000
6.9873238398 49.9280496142 -3.0462 0.0000
Testing
As a test, I then went back to see how far each IP was from what I thought its ideal location was.
The result is that all but three of the 13 points are within 20 metres.
(Plat C-46 is unusually far off, at 48 metres. I took a second look at it. Privately, I wonder if they didn’t mistakenly exchange the final digits of the easting and northing coordinates on this one. The coordinates given on the plat are E 361.56, N 239.43, but E 361.53, N 239.46 would have put it right where I expected it. But that’s just me trying to massage the data to fit my hypothesis!)
In the end, DoughboyMIA is using the cct command-line given above, and it seems to be going well. Robert Laplander, who heads up the project of teams working on the ground in France, says that they used/tested the process in the summer of 2023, “and it worked great.”
I don’t know whether the Helmert parameters could be further refined. It would be nice to have a much bigger dataset of Lambert coordinates and their correct corresponding WGS84 coordinates, and see the mean solution which that provides. But there’s also the consideration that those military surveyors must have sometimes made mistakes, and about a quarter of the time the Lamberts are rounded to the nearest 100 metres. So we may not be able to get much better than this.
In addition to the history, this process taught me a huge amount about coordinate transformations. I used to think we knew how to how to transform between any two projections, but we don’t. Now when QGIS tells me it used a “ballpark” transformation, warning bells go off!
The Wandering Cartographer 