Под катом будет небольшая заметка о применении пространственного индекса
на основе zcurve для индексации точечных данных, расположенных на сфере. А так же bencmark-и для PostgreSQL и сравнение с таким же (но совсем другим)
индексом на R-дереве.

В развитие темы (1, 2, 3, 4, 5, 6).
Собственно, возвращаемся к самому началу — идее индексировать географические координаты, размещая их на поверхности сферы. Обычная индексация пары широта/долгота приводит к искажениям вдали от экватора, работа с проекциями не универсальна. Поэтому мысль переводить географические координаты в трехмерное пространство выглядит довольно изящно.

Сама по себе идея эта не нова, аналогично работает, например, расширение PostgreSQL — PGSphere, которое использует для индексации 3-мерное R-дерево. С ним и будем сравнивать.

Подготовка данных.

PGSphere

• Для начала придётся выкачать, собрать и инсталлировать расширение (автор использовал текущую версию 1.1.5)
  
  

  gmake USE_PGXS=1 PG_CONFIG=/usr/bin/pg_config
  sudo gmake USE_PGXS=1 PG_CONFIG=/usr/bin/pg_config install
  

• Далее загрузить его (psql)
  
  

  CREATE EXTENSION pg_sphere;

• Создадим таблицу для тестовых данных
  
  

  CREATE TABLE spoint_data (sp spoint);

• Нам потребуется источник случайных данных.
  Первый параметр программы — радиус, второй — число результатов.
  Единственная тонкость — данные равномерно распределены внутри шара с заданным радиусом, иначе не получится равномерного распределения на сфере
• Случайные данные пропустим через скрипт awk чтобы превратить в геокоординаты
  
  

  # --- gendata.awk ------
  BEGIN{
  	pi=3.1415926535897932;
  	degra=pi/180.0;
  	rad=180.0/pi;
  	Grad = 1000000.;
  }
  {
  	x = $1;	y = $2;	z = $3;
  	r3 = sqrt(x*x + y*y + z*z);
  	x *= Grad / r3;
  	y *= Grad / r3;
  	z *= Grad / r3;
  
  	r2 = sqrt(x*x + y*y);
  	lat = atan2(z, r2) * rad;
  	lon = 180. + atan2(y, x) * rad;
  
  	printf ("(%14.10fd, %14.10fd)\n", lon, lat);
  }

• Собственно создание данных, здесь радиус не важен, важно чтобы и pgsphere и zcurve получили одни и те же данные. Сортировка весьма желательна для ускорения индексации.
  
  

  ./random 1000000 100000000 | gawk -f gendata.awk | sort > pgsphere.txt

• Заливаем данные в нашу таблицу
  
  

  COPY spoint_data (sp) FROM  '/home/.../pgsphere.txt';

• Индексируем
  
  

  CREATE INDEX sp_idx ON spoint_data USING gist (sp);

ZORDER

• Для начала придётся выкачать, собрать и инсталлировать расширение
  
  

  gmake USE_PGXS=1 PG_CONFIG=/usr/bin/pg_config
  sudo gmake USE_PGXS=1 PG_CONFIG=/usr/bin/pg_config install
  

• Создадим таблицу для тестовых данных
  
  

  create table test_points_3d (x integer,y integer, z integer);

• Нам потребуется тот же источник случайных данных.
• Случайные данные пропустим через скрипт awk чтобы разместить их внутри куба со стороной в 2 000 000
  
  

  #--- gendata2.awk ------
  BEGIN{
  	pi=3.1415926535897932;
  	degra=pi/180.0;
  	rad=180.0/pi;
  	Grad = 1000000.;
  }
  {
  	x = $1;	y = $2;	z = $3;
  	r3 = sqrt(x*x + y*y + z*z);
  	x *= Grad / r3;
  	y *= Grad / r3;
  	z *= Grad / r3;
  
  	ix = int(x+0.5+Grad);
  	iy = int(y+0.5+Grad);
  	iz = int(z+0.5+Grad);
  	print ix"\t"iy"\t"iz;
  }

• Собственно создание данных, здесь радиус важен. Сортировка не обязательна.
  
  

  ./random 1000000 100000000 | gawk -f gendata2.awk > zcurve.txt

• Заливаем данные в нашу таблицу
  
  

  COPY test_points_3d FROM '/home/.../zcurve.txt';

• Индексируем
  
  

  create index zcurve_test_points_3d on test_points_3d(zcurve_num_from_xyz(x,y,z));

Подготовка тестов

PGSphere

Для тестирования потребуется вот такой awk скрипт

#--- gentest.awk -------
BEGIN{
	pi=3.1415926535897932;
	degra=pi/180.0;
	rad=180.0/pi;
	Grad = 1000000.;
}
{
	x = $1;	y = $2;	z = $3;
	r3 = sqrt(x*x + y*y + z*z);
	x *= Grad / r3;
	y *= Grad / r3;
	z *= Grad / r3;

	r2 = sqrt(x*x + y*y);

	lat = atan2(z, r2) * rad;
	lon = 180. + atan2(y, x) * rad;

#	EXPLAIN (ANALYZE,BUFFERS) 
  printf ("select count(1) from spoint_data where sp        @'<(%14.10fd,%14.10fd),.316d>'::scircle;\n", lon, lat);
}

Этот скрипт вполне симметричен тому, с помощью которого мы подготавливали данные. Стоит обратить внимание на число .316, это радиус сферы с центром в вычисленной случайной точке, в которой мы ищем данные

Подготовка серии запросов делается так:

./random 1000000 100 1023 | gawk -f gentest.awk >tests1.sql

Здесь 100 — размер тестовой серии, 1023 — seed рандомизатора.

ZCURVE

Для тестирования тоже потребуется awk скрипт

#--- gentest2.awk -------
BEGIN{
	pi=3.1415926535897932;
	degra=pi/180.0;
	rad=180.0/pi;
	Grad = 1000000.;
}
{
	x = $1;	y = $2;	z = $3;
	r3 = sqrt(x*x + y*y + z*z);
	x *= Grad / r3;
	y *= Grad / r3;
	z *= Grad / r3;

	ix = int(x+0.5+Grad);
	iy = int(y+0.5+Grad);
	iz = int(z+0.5+Grad);
#	EXPLAIN (ANALYZE,BUFFERS) 
	lrad = int(0.5 + Grad * sin(.316 * degra));
	print "select count(1) from zcurve_3d_lookup_tidonly('zcurve_test_points_3d',   "ix-lrad","iy-lrad","iz-lrad","ix+lrad","iy+lrad","iz+lrad");";
}

Этот скрипт вполне симметричен тому, с помощью которого мы подготавливали данные. Опять же, обращаем внимание на число .316, это половина стороны куба с центром в вычисленной случайной точке, в котором мы ищем данные.

Подготовка серии запросов делается так:

./random 1000000 100 1023 | gawk -f gentest2.awk >tests2.sql

Здесь 100 — размер тестовой серии, 1023 — seed рандомизатора, лучше, если он совпадает с оным от pgsphere.

Benchmark

Как и раньше, замеры проводились на скромной виртуальной машине с двумя ядрами и 4 Гб ОЗУ, поэтому времена не имеют абсолютной ценности, а вот числам прочитанных страниц по прежнему можно доверять.

Времена показаны на вторых прогонах, на разогретом сервере и виртуальной машине. Количества прочитанных буферов — на свеже-поднятом сервере.

Radius	AVG NPoints	Nreq	Type	Time(ms)	Reads	Hits
.01°	1.17 0.7631 (0.7615)	10 000	zcurve rtree	.37 .46	1.4397 2.1165	9.5647 3.087
.0316°	11.6 7.6392 (7.6045)	10 000	zcurve rtree	.39 .67	2.0466 3.0944	20.9707 2.7769
.1°	115.22 76.193 (76.15)	1 000	zcurve rtree	.44 2.75 *	4.4184 6.073	82.8572 2.469
.316°	1145.3 758.37 (760.45)	1 000	zcurve rtree	.59 18.3 *	15.2719 21.706	401.791 1.62
1.°	11310 7602 (7615)	100	zcurve rtree	7.2 94.5 *	74.9544 132.15	1651.45 1.12

где
    Radius — размер поисковой области в градусах
    Npoints — среднее число точек в выдаче, в скобках — теоретически ожидаемое число
     (в сфере 41252.96 кв. градусов, 100 000 000 точек, ~2424 точки на кв. градус)
    Nreq — число запросов в серии
    Type —
      ‘zcurve’ — оно и есть
      ’rtree’- PGSphere
    Time(ms) — среднее время выполнения запроса
    Reads — среднее число чтений на запрос
    Hits — число обращений к буферам

    * в какой-то момент производительность R-tree начинает резко
    проседать, связано это с тем, это дерево читает заметно больше
    страниц и его рабочий набор перестаёт помещаться в кэше (по-видимому).

Отметим, что zcurve находит больше данных, что логично т.к. он ищет внутри куба, а не сферы как PGSphere. Поэтому требуется пост-фильтрация в духе HAVERSINE. Но здесь мы сравниваем только производительности индексов.

Попробуем оценить разницу. Вообще, задача найти площадь куска сферы, попавшей внутрь куба нетривиальна. Попробуем сделать оценку.

• Предположим, что наш куб в среднем вырезает из сферы ту же площадь, что и сфера равного объема
• Объем единичной сферы 1.33*3.14=4.19
• Объем куба со стороной 2 = 8.
• Тогда корень третьей степени из 8/4.19 = 1.24 — это отношение радиусов мнимой сферы к настоящей
• соотношение площадей мнимой сферы к настоящей 1.24*1.24=1.54
• имеем из экспериментальных данных 1.17/0.7631= 1.5332
  Bingo!