2022 / Day 19: Globe
Back almost 15-20 years ago when I started working in GIS I remember an older colleague telling me that when choosing coordinate precision in degrees for the area at the latitudes of Estonia the rule of thumb is approximately:
- 1° for 100 km
- 0.1° for 10 km
- 0.01° for 1 km / 1000 m
- 0.001° for 100 m
- 0.0001° for 10 m
- 0.00001° for 1 m
- 0.000001° for 0.1 m / 10 cm
- 0.0000001° for 1 cm
- and so on
Not exact figures, but more or less.
This is exactly what came to my mind when some time ago the xkcd coordinate precision comic was making its rounds
Coordinate Precision https://t.co/6YMvsNzpDa https://t.co/xZssjYfNGZ pic.twitter.com/UJfTnMuycS
— XKCD Comic (@xkcdComic) July 1, 2019
The funny (at least to me in my mind) thing with a reference system in degrees on a sphere is that the further away you move from the equator towards any of the poles the less precision you would need for your data not to be packed into a “grid”. So for example (but depending on what you’re trying to do - see the xkcd comic above) on the equator it makes a difference if you use 4 decimal places or 5 decimal places for coordinates in degrees. But if you’re exactly at the pole it makes no difference whatsoever - when the latitude is 90° (or -90°) it doesn’t really matter if the precision is zero, one, five, or whatever decimal places. Even more - the meridian value will not have any impact at all. You’ll be -180° .. 180° west and east all the time and at the same time.
Today’s entry to the challenge depicts the relationship between ground distance unit (meter) and the longitude coordinate decimal degree precision as we move from the equator towards any of the poles. Latitude lines of 20 of the worlds biggest cities (according to Wikipedias List of largest cities with their locations from NaturalEarth) are also shown. And in addition Longyearbyen and Ushuaia as a reference. These will serve as the y-axis grid.
The isoprecision lines themselves are calculated in the CTE pts by creating
a line from the Null Island to the North pole along the Greenwich meridian
with st_makeline
and adding vertices to it every 0.1 degrees with
st_segmentize and dumping
all vertices with st_dumppoints.
These will serve as the “original points”.
Next st_translate all the points
by delta longitude from a series of precisions in translates.
And then create the isoprecision line by creating new nodes at
the “original point” st_y
but setting the x-coordinate at the value of
st_distancespheroid
between the “original point” and the corresponding translated point. Finally
aggregate these together into a linestring over the used precision value and
in the order of original points.
with
pts as (
select
(
st_dumppoints(
st_segmentize(
st_makeline(
st_point(0,0,4326),
st_point(0,90,4326)
),
0.1
)
)
).*
),
translates as (
select
s, pts.path[1] as ord,
pts.geom as original, st_translate(pts.geom, s, 0) as geom
from
pts,
unnest(
array[
0.001,
0.0001,
0.00001,
0.000001
]
) s
),
cities as (
select
name, cast(lon as numeric) lon,
cast(lat as numeric),
round(cast(lat as numeric),1) as rounded_lat
from (
values (
'Longyearbyen', 15.63, 78.217
),(
'Ushuaia', -68.3, -54.8
),(
'North/South pole', 0.0, 90.0
),(
'Equator', 0.0, 0.0
),(
'Tokyo', 139.749462, 35.686963
),(
'Delhi', 77.199980, 28.600023
),(
'Shanghai', 121.434559, 31.218398
),(
'São Paulo', -46.626966, -23.556734
),(
'Mexico City', -99.132934, 19.444388
),(
'Cairo', 31.248022, 30.051906
),(
'Mumbai', 72.855043, 19.018936
),(
'Beijing', 116.386340, 39.930838
),(
'Dhaka', 90.406634, 23.725006
),(
'Osaka', 135.458199, 34.751981
),(
'New York', -73.981963, 40.751925
),(
'Karachi', 67.01, 24.86
),(
'Buenos Aires', -58.399477, -34.600556
),(
'Chongqing', 106.5504, 29.5637
),(
'Istanbul', 29.008056, 41.106942
),(
'Kolkata',88.37, 22.5675
),(
'Manila',120.980271, 14.606105
),(
'Lagos',3.389585,6.445208
),(
'Rio do Janeiro',-43.226967, -22.923077
),(
'Tianjin', 117.2054, 39.1336
)
) x (name, lon, lat)
),
lines as (
select
row_number() over()::int as oid,
s::varchar,
st_makeline(
array_agg(
st_point(
st_distancespheroid(
original,
geom,
'SPHEROID["WGS 84",6378137,298.257223563]'
),
st_y(original)
) order by ord
)
) as geom
from translates f
group by s
)
select
row_number() over():: int as id, oid, name, geom, cl
from (
select
oid, name,
case
when mod(oid,2)=0 then
st_makeline(
st_point(-2.0, abs(c.rounded_lat)),
st_point(130.0, abs(c.rounded_lat))
)
else
st_makeline(
st_point(-10.0, abs(c.rounded_lat)),
st_point(122.0, abs(c.rounded_lat))
)
end as geom,
'cities' as cl
from (
select
row_number() over(order by abs(rounded_lat))::int as oid,
c.name, c.rounded_lat
from
cities c
) c
union all
select
oid, s as name, geom, 'isoprecision' as cl
from
lines
union all
select
null, m::varchar as name,
st_makeline(st_point(m,-5), st_point(m,92)),
'meter-grid' as cl
from
generate_series(0, 120, 5) m
) d
;