DEV Community

Srdjan Popovic
Srdjan Popovic

Posted on

Same Zone, Same Projection, 427 Metres Apart

#gis #postgres #geospatial #datascience

I was writing a blog post about the architecture of a system I built when I found a bug in it. This is that bug, and the much more interesting thing standing behind it.

The system ingests mobile mapping surveys — LiDAR trajectories with spherical photos along them, driven along roads. The files arrive in whatever coordinate system the surveyor was working in, which in this part of the world is not a settled question. So the upload form asks the user to pick one:

SOURCE_SRID_CHOICES = [
    (32634, 'UTM zone 34N (EPSG:32634)'),
    (31277, 'Gauss-Kruger zone 7 (EPSG:31277)'),
    (8686,  'Gauss-Kruger zone 7 (EPSG:8686)'),
]
Enter fullscreen mode Exit fullscreen mode

Two codes, one label. That was deliberate — I knew there were two definitions floating around for the same zone and wanted users to be able to pick the one matching their source.

I was wrong about what the second one was.

What the database says

You don't have to guess about an EPSG code. If you're running PostGIS, the definition is sitting in spatial_ref_sys and you can just read it.

SELECT srid, proj4text FROM spatial_ref_sys WHERE srid IN (31277, 8686);
Enter fullscreen mode Exit fullscreen mode 
31277 | +proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0
        +datum=hermannskogel +units=m +no_defs

 8686 | +proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=500000 +y_0=0
        +ellps=bessel +towgs84=476.08,125.947,417.81,... +units=m +no_defs
Enter fullscreen mode Exit fullscreen mode

Look at lon_0. One is on the 21st meridian, the other on the 15th. Look at x_0: 7,500,000 versus 500,000.

These are not two definitions of the same zone. 8686 is MGI 1901 / Slovenia Grid — a different country's national grid, six degrees west.

How far wrong? Take a coordinate pair from Belgrade and push it through both:

SELECT ST_AsText(ST_Transform(ST_SetSRID(ST_MakePoint(7457000, 4958000), 31277), 4326));
-- POINT(20.456 44.765)   Belgrade. Correct.

SELECT ST_AsText(ST_Transform(ST_SetSRID(ST_MakePoint(7457000, 4958000), 8686), 4326));
-- POINT(76.635 25.100)   Rajasthan, India. About 5,400 km off.
Enter fullscreen mode Exit fullscreen mode

Nobody had used the option yet — every record in the table had the field null — so this cost nothing except my confidence. But it's in a dropdown, labelled invitingly, waiting.

Why that mistake was easy to make

Here's the part that turns a careless bug into something worth writing about.

The modern EPSG codes for this family are the MGI 1901 / Balkans zone N series. If you go looking for them in spatial_ref_sys, you get this:

SELECT srid, srtext FROM spatial_ref_sys WHERE srtext LIKE '%MGI 1901%Balkans%';
Enter fullscreen mode Exit fullscreen mode 
3907  MGI 1901 / Balkans zone 5   lon_0=15
3908  MGI 1901 / Balkans zone 6   lon_0=18
3909  MGI 1901 / Balkans zone 7   lon_0=21
3910  MGI 1901 / Balkans zone 8   lon_0=24

8677  MGI 1901 / Balkans zone 5   lon_0=15
8678  MGI 1901 / Balkans zone 6   lon_0=18
8679  MGI 1901 / Balkans zone 8   lon_0=24
Enter fullscreen mode Exit fullscreen mode

Read the second block again. 8677 is zone 5, 8678 is zone 6, and 8679 is zone 8.

The 867x series skips zone 7. There is no 8679 = zone 7 to find, because zone 7 already had 3909 and the later batch didn't re-issue it. So anyone scanning that numeric range for "the modern zone 7 code" finds a gap where their pattern-matching brain expects a hit, and the nearby numbers are all real, plausible, MGI-1901-family codes.

8686 is one of those nearby numbers. It is a real code. It belongs to the same datum family. It is in the same region of Europe. It is simply the wrong country.

The correct code for Gauss-Krüger zone 7 in this region is EPSG:3909.

The trap that actually matters

Fixing the label is a one-line change. The thing worth your attention is the other pair in that dropdown, because it fails in a much more dangerous way.

31277 and 3909 both describe Gauss-Krüger zone 7. Same central meridian, same scale factor, same false easting, same ellipsoid:

3909   +proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +ellps=bessel
       +towgs84=682,-203,480,0,0,0,0

31277  +proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000
       +datum=hermannskogel
Enter fullscreen mode Exit fullscreen mode

Every projection parameter is identical. The only difference is how the datum reaches WGS 84: an explicit three-parameter shift in one, PROJ's built-in Hermannskogel datum in the other.

Same input through both:

3909   POINT(20.451342 44.765252)
31277  POINT(20.456732 44.765230)
Enter fullscreen mode Exit fullscreen mode 
-- distance between them
426.73 metres
Enter fullscreen mode Exit fullscreen mode

427 metres. From an identical-looking projection definition, differing only in a datum parameter you have to scroll sideways to notice.

And here is why that is worse than the 5,400 km error: it looks fine.

Five thousand kilometres announces itself. You load the layer, it isn't on the map, you find the problem in thirty seconds. Four hundred metres puts your data in the right country, the right city, the right neighbourhood. Zoom out to a national view and both versions are the same pixel. Everything renders. Nothing errors.

It only becomes visible when someone zooms in far enough to notice that the road inventory is on the wrong side of the road — in a dataset whose entire value proposition is sub-metre positional accuracy.

Why there are two in the first place

Briefly, because it explains why you can't just delete the old one.

31277 is marked deprecated in the EPSG registry. It's the older realization, and PROJ resolves its datum through Hermannskogel. 3909 is the current code with an explicit shift to WGS 84.

But deprecated doesn't mean unused. Decades of surveying in this region was done, stored, and delivered in the older definition, and files still arrive that way. If you only offer the modern code, you silently apply the wrong datum shift to legacy data — the same 427 m error, in the other direction.

So both belong in the dropdown. What doesn't belong is giving them the same label.

What I changed

The codes stay. The labels stop lying:

SOURCE_SRID_CHOICES = [
    (32634, 'UTM zone 34N — WGS 84 (EPSG:32634)'),
    (3909,  'Gauss-Krüger zone 7 — MGI 1901 (EPSG:3909)'),
    (31277, 'Gauss-Krüger zone 7 — MGI/Hermannskogel, legacy (EPSG:31277)'),
]
Enter fullscreen mode Exit fullscreen mode

8686 is gone. 3909 replaces it. The two remaining zone 7 entries name their datum, because the datum is the entire difference between them.

What I'd take from this

Read spatial_ref_sys instead of trusting the label. The definition is in your database. SELECT proj4text takes five seconds and tells you the central meridian, the false easting, and the datum. Every mistake in this post was visible in that one column.

Compare lon_0 and x_0 first. They're the parameters that produce catastrophic, obvious errors, so they're the cheapest to check. If those match and you still have two codes, your difference is in the datum — which is the expensive, quiet kind.

Add a plausibility check on transformed coordinates. One bounding box test — does the result land inside the country this data is supposed to be in — would have caught the Slovenia Grid mistake automatically, at upload time, before it ever reached a map. It would not have caught the 427 m one, which is the point: cheap checks catch the loud failures, and you need to know they don't catch the quiet ones.

Never label two coordinate systems identically. If the user has to choose between them, the label has to contain the thing that differs. "Gauss-Krüger zone 7" twice is not a choice, it's a coin flip with a 427 m stake.

Gap-filled numbering, deprecated-but-still-in-use codes, and neighbouring countries sharing a datum family — none of these are anyone's fault. They're the residue of a century of national surveying, encoded into a flat integer namespace. The registry is doing its job. The label was doing mine, badly.

Top comments (0)