pycairo и сплайны

24 июня, 2009

Вот тут есть неплохой пост про использование pycairo для сглаживания ломанных линий в сплайны Безье. Вот только код не работает если попытаться сгладить вот такую ломанную:ломанная

curve_

Правильно нарисован только верхний ряд. Нижний и вертикальный сгладились как-то шиворот-навыворот.

Исправил функцию prepare_curve_data с использованием библиотеки vec2d

def prepare_curve_data(graph_data):
    from vec2d import *
    prepared_data = []
    for i in range(len(graph_data)):
        x, y = graph_data[i][0], graph_data[i][1]
        cur_point = vec2d(graph_data[i][0], graph_data[i][1])
        
        
        if (i != 0) and (i != len(graph_data) - 1):
            back_point = vec2d(graph_data[i - 1][0], graph_data[i - 1][1])
            forw_point = vec2d(graph_data[i + 1][0], graph_data[i + 1][1])
                        
            back_vect = back_point - cur_point
            forw_vect = forw_point - cur_point
            back_vect_proj = back_vect / 2.0
            forw_vect_proj = forw_vect / 2.0
            
            angle_between = back_vect.get_angle_between(forw_vect) / 2.0          
            
            back_vect.rotate(angle_between - 90)
            forw_vect.rotate(90 - angle_between)
            
            back_vect = back_vect_proj.projection(back_vect)
            forw_vect = forw_vect_proj.projection(forw_vect)
            
            back_cpoint = cur_point + back_vect
            forw_cpoint = cur_point + forw_vect
            
            cx1, cy1 = back_cpoint[0], back_cpoint[1]
            cx2, cy2 = forw_cpoint[0], forw_cpoint[1]
        
        else:
           cx1, cy1, cx2, cy2 = x, y, x, y 
        
        prepared_data.append((x, y, cx1, cy1, cx2, cy2))
        
    return prepared_data

И вот результат:

curve

Программа, рисующая эту картинку:


import cairo
from math import pi, sqrt

width = 500
height = 500
graph_data = [
(0, 10),(20, 50),(40, 80),(60, 5),(80, 10),(100, 20),
(120, 30),(140, 60),(160, 95),(180, 30),(200, 50),
(220, 70),(240, 80),(260, 10),(280, 60),(300, 30),
(320, 90),(340, 95),(360, 30),(380, 10),(400, 5),
(420, 20),(440, 80),(460, 70),(480, 20),(500, 40),
(490, 0), (450, 20), (420, 40), (495, 60), (490, 80),
(480, 100), (470, 120), (440, 140), (405, 160), 
(470, 180), (450, 200), (430, 220), (420, 240),
(490, 260), (440, 280), (470, 300), (410, 320),
(405, 340), (470, 360), (490, 380), (495, 400), 
(480, 420), (420, 440), (430, 460), (480, 480), (460, 500),
(500, 490), (480, 450), (460, 420), (440, 495), (420, 490), 
(400, 480), (380, 470), (360, 440), (340, 405), (320, 470), 
(300, 450), (280, 430), (260, 420), (240, 490), (220, 440), 
(200, 470), (180, 410), (160, 405), (140, 470), (120, 490), 
(100, 495), (80, 480), (60, 420), (40, 430), (20, 480), (0, 460)
]

def prepare_curve_data(graph_data):
    from vec2d import *
    prepared_data = []
    for i in range(len(graph_data)):
        x, y = graph_data[i][0], graph_data[i][1]
        cur_point = vec2d(graph_data[i][0], graph_data[i][1])
        
        
        if (i != 0) and (i != len(graph_data) - 1):
            back_point = vec2d(graph_data[i - 1][0], graph_data[i - 1][1])
            forw_point = vec2d(graph_data[i + 1][0], graph_data[i + 1][1])
                        
            back_vect = back_point - cur_point
            forw_vect = forw_point - cur_point
            back_vect_proj = back_vect / 2.0
            forw_vect_proj = forw_vect / 2.0
            
            angle_between = back_vect.get_angle_between(forw_vect) / 2.0          
            
            back_vect.rotate(angle_between - 90)
            forw_vect.rotate(90 - angle_between)
            
            back_vect = back_vect_proj.projection(back_vect)
            forw_vect = forw_vect_proj.projection(forw_vect)
            
            back_cpoint = cur_point + back_vect
            forw_cpoint = cur_point + forw_vect
            
            cx1, cy1 = back_cpoint[0], back_cpoint[1]
            cx2, cy2 = forw_cpoint[0], forw_cpoint[1]
        
        else:
           cx1, cy1, cx2, cy2 = x, y, x, y 
        
        prepared_data.append((x, y, cx1, cy1, cx2, cy2))
        
    return prepared_data

def draw_point(x, y, opacity):
    cr.move_to(x + 2, y)
    cr.arc(x, y, 2, 0, 2 * pi)
    cr.set_source_rgba(0, 0, 1, opacity)
    cr.stroke()

def debug_points(cr, prepared_data):
    for i in range(0, len(prepared_data)):
        x, y = prepared_data[i][0], prepared_data[i][1]
        cx1, cy1 = prepared_data[i][2], prepared_data[i][3]
        cx2, cy2 = prepared_data[i][4], prepared_data[i][5]

        draw_point(x, y, 1)
        if cx1 != x or cy1 != y:
            draw_point(cx1, cy1, 0.3)
            cr.move_to(x, y)
            cr.line_to(cx1, cy1)
            cr.set_source_rgb(0.5, 0.5, 0.5)
            cr.stroke()

        if cx2 != x or cy2 != y:
            draw_point(cx2, cy2, 0.3)
            cr.move_to(x, y)
            cr.line_to(cx2, cy2)
            cr.set_source_rgb(0.5, 0.5, 0.5)
            cr.stroke()

def poly_curve(cr, prepared_data):
    for i in range(0, len(prepared_data) - 1):
        x, y = prepared_data[i][0], prepared_data[i][1]
        cx1, cy1 = prepared_data[i][4], prepared_data[i][5]
        cx2, cy2 = prepared_data[i + 1][2], prepared_data[i + 1][3]
        x2, y2 = prepared_data[i + 1][0], prepared_data[i + 1][1]
        cr.move_to(x, y)
        cr.curve_to(cx1, cy1, cx2, cy2, x2, y2)


surface = cairo.ImageSurface(cairo.FORMAT_ARGB32, width, height)
cr = cairo.Context(surface)
cr.set_line_width(1)
cr.set_source_rgb(1, 1, 1)
cr.set_operator (cairo.OPERATOR_SOURCE)
cr.paint()

prepared_data = prepare_curve_data(graph_data)
debug_points(cr, prepared_data)
cr.set_line_width(2)
poly_curve(cr, prepared_data)
cr.set_source_rgb(0, 0, 0)
cr.stroke()

surface.write_to_png('curve.png')


Построение изолиний в python

28 мая, 2009

Понадобилось мне для матрицы значений получить координаты изолинии этих значений. Матрица, разумеется, с регулярным шагом. Так как основная моя программа написана на python, то и программу, реализующую это тоже хотелось бы на питоне. Однако стихийный поиск не навёл меня на уже готовое решение. Разумеется, это можно было бы реализовать при помощи PyNGL или тому подобного, но мне необходимы сами координаты этих изолиний.
В общем, после непродолжительных мытарств, решился написать свой велосипед. Задача решена в лоб, и примитивно. Но, может, кому пригодиться.

Пример использования:

import isolines

data = [[21, 35, 55, 41], \
[24, 47, 77, 56], \
[12, 46, 54, 27], \
[8, 70, 82, 38], \
[6, 81, 109, 63], \
[9, 56, 77, 47], \
[11, 27, 54, 46], \
[20, 38, 82, 70]]

# список интервалов, через которые будут проведены изолинии
intervals=range(0,100,10)

»’
prepared_matrix — подготовленные ячейки, для дальнейших расчетов
calculated_dots — точки пересечения граней квадратов с горизонтальной плоскостью
coords — результат. координаты отрезков изолиний.
»’
prepared_matrix=isolines.prepare_matrix(data)
calculated_dots=isolines.calculate_dots(prepared_matrix,intervals)
coords=isolines.calculate_coords(calculated_dots)

print coords

Результат выполнения:

[[[30, 1.1428571428571428, 0.5, 0.76086956521739135, 1.5], \
[40, 1.5, 0.91666666666666674, 1.1956521739130435, 1.5]], \
[[40, 1.75, 0.5, 1.5, 0.91666666666666674], \
[50, 2.25, 0.5, 1.6000000000000001, 1.5], \
[60, 2.5, 0.72727272727272729, 1.9333333333333333, 1.5], \
[70, 2.5, 1.1818181818181817, 2.2666666666666666, 1.5]], ….

Список, в котором отдельно для каждой ячейки и для каждой стороны этой ячейки содержатся координаты отрезков изолиний, и их значений в виде:
[значение изолинии, x1, y1, x2, y2]

А вот как это решение выглядит в графическом виде:
Изолинии
Не самая красивая диаграмма. Но характер распределения передает.

Может, у кого-то есть более изящное решение?


ВНИМАНИЕ!

1 апреля, 2009

Недавно исследователи открыли факт заражения наших водопроводных систем опасным химикатом. Этот химикат бесцветный, безвкусный и не имеет запаха. Правительство не предприняло никаких попыток регулирования этого опасного заражения. Данный химикат называется дигидрогена монооксид (Dihydrogen monoxide).

Химикат используется для следующих целей:
В производстве как растворитель и охладитель
В ядерных реакторах
В производстве пенопласта
В огнетушителях
В химических и биологических лабораториях
В производстве пестицидов
В искусственных пищевых добавках
Химикат является основной составляющей кислотных дождей
Действует на эрозию почвы
Ускоряет коррозию и вредит большинству электроприборов
Длительный контакт с химикатом в его твёрдой форме приводит к серьёзным повреждениям кожи человека
Контакт с газообразной формой химиката приводит к сильным ожогам
Вдыхание даже небольшого количества химиката грозит смертельным исходом
Химикат обнаружен в злокачественных опухолях, нарывах, язвах и прочих болезненных изменениях тела
Химикат развивает наркозависимость; жертвам при воздержании от потребления химиката грозит смерть в течение 168 часов
Ни один известный очиститель не способен полностью очистить воду от этого химиката.

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

Мы призываем население проявить сознательность и протестовать против дальнейшего использования этого опасного химиката.

Источник http://ru.wikipedia.org/wiki/Дигидроген_моноксид


Построение стереограмм в linux

12 августа, 2008

Построение стереограмм трещиноватости является повседневной задачей в структурной геологии. Для её решения уже разработано большое количество программных средств под Windows, но под Linux специальной утилиты пока не существует. Есть лишь небольшое количество малофункциональных программ.

Но так как строить такие стереограммы мне необходимо, да еще и в больших количествах, я написал скрипт, осуществляющий такие построения, используя функционал GMT.

#!/bin/bash
fig=test.ps
title=»Проверка»
grid=test.grd
cpt=test.cpt

gmtset CHAR_ENCODING ISO-8859-5

echo ‘0, -90’ | psxy -R0/360/-90/0 \
-JS0/-90/16c -Sx0.4c -Ba90f90/g90:.$title: \
—BASEMAP_TYPE=plain -K -V > $fig

#рисуем нормали к плоскостям
cat test_n.txt | awk -F»\t» ‘{print $1, $2-90}’ \
| psxy -R -J -Sc0.2c -N -Gblack -Wwhite -O -K >> $fig

#рисуем направления смещений по плоскостям
cat test_v.txt | awk -F»\t» ‘{print $1, $2-90, -$3+90, «1»}’ \
| psxy -R -J -SV0.08/0.3c/0.12c -N -Gblack -O -K >> $fig

echo «50 5 12 0 0 CT Тестовая стереограмма» | pstext -R0/50/0/50 -JX16c -N -O >> $fig

Результатом выполнения скрипта будет вот такая стереограмма:


Пространственная кластеризация точечных данных

12 мая, 2008

Пост является переводом Spatial Clustering of Point Data с небольшими дополнениями.

Задача: Есть n точек, распределенных каким-то образом в пределах заданной территории. Каждая точка кроме пространственного положения обладает значением (например, содержанием химического элемента)
Требуется объединить точки в m групп.

Инструменты: GRASS, R.

Задача теоретическая. Точки генерируем случайным образом, и по ходу решения задачи, попытаемся определиться с числом групп в которые эти точки надо объединить.
Для наглядности сгенерируем растровую карту и возьмем с неё значения для точек:

r.surf.fractal out=surf.fract d=2.05
v.random out=pnts n=100
v.drape pnts type=point rast=surf.fract out=points

d.rast surf.fract
d.vect points icon=basic/box fcol=black col=black size=6

v.out.ascii in=points out=points.xy

На соседнем терминале запускаем R

x <- read.table(‘points.xy’,’sep=’|’)
names(x) <- c(‘easting’, ‘northing’, ‘surf’ ,’cat’)
# отсекаем столбец cat
y <- data.frame(x[,1:3])
row.names(y) <- x$cat

# загружаем библиотеки для кластерного анализа
library(cluster)
library(flexclust)

# попытаемся найти оптимальное число групп
s <- stepFlexclust(y, K=2:10, nrep=20)
plot(s)

# похоже, оптимальное число 5
y.pam <- pam(y, 5, stand=TRUE)

# подготовим данные для экспорта
y$cluster <- y.pam$clustering
y$orig_cat <- as.numeric(row.names(y))

# экспортируем данные в текстовый файл
write.table(y, file=’points.clust’, row.names=FALSE)

возвращаемся в GRASS и импортируем текстовый файл

v.in.ascii in=points.clust out=pclust fs=» » columns=’x double, y double, srf double, orig_cat integer, cluster integer’ skip=1

for x in $(seq 1 5)
do v.extract —o in=pclust where=»cluster=$x» out=pclust_$x
v.hull —o in=pclust_$x out=pclust_hull_$x
d.vect pclust_hull_$x type=boundary fcol=none width=2 col=white
done
d.vect pclust icon=basic/box fcol=black col=black size=6


Начало работы с GMT

19 марта, 2008

На сайте gis-lab опубликована статья «Начало работы с GMT«, в которой показаны основы работы с этим пакетом. Спасибо Максиму Дубинину за помощь и терпеливую корректировку сего опуса.


Русское руководство по GRASS

4 марта, 2008

Наконец-то, после напряженных месяцев кропотливой работы, ребята с сайта gis-lab.info завершили перевод на русский язык учебного руководства по GRASS. Молодцы, большое дело сделали!


Почему работа в командной строке — умирающее искусство?

6 февраля, 2008

На днях на блоге PerryGeo обнаружил сей замечательный пост. Думаю, он актуален и для многих отечественных ГИС (и не только) пользователей. Поэтому, ниже предлагаю вольный перевод:

К сожалению, большинство пользователей ГИС довольно далеки от командной строки (CLI — command line interface). Порой, даже опытные пользователи, оказываясь один на один с командной строкой, испытывают шок. Насквозь оконные интерфейсы нынешних лидеров ГИС-рынка, облекающие любую операцию в ГУЕвую оболочку (GUI — graphical user interface), только способствуют этому (тех кто еще помнит как работать в командной строке ESRI Arc/Info уже называют «старая гвардия»). К тому же, у пользователей XP и Vista отсутствует доступ к командной строке DOS. У Linux-пользователей дела с этим обстоят по-лучше, но разница уже не так велика, по причине наступления таких дистрибутивов, как Ubuntu (с через чур ‘дружелюбным’ интерфейсом).
Так что же такого ужасного в командной строке? Почему считается, что командная строка сложнее графического интерфейса? Я пришел к выводу что в некоторых случаях все совсем наоборот — набрать что-то и получить назад ответ проще простого! К тому же, создается ощущение полного контроля за компьютером (что, на самом деле правда). Компьютер всегда выполняет только те приказания, которые вы ему даете, каким интерфейсом вы бы не пользовались (GUI или CLI). Вот только GUI не концентрируется на мелочах, поэтому, вам не обязательно точно знать что вы приказали компьютеру. Такая легкость дается ценой многих важных факторов. Читать далее…


Установка GMT в Debian/Ububtu

1 февраля, 2008

Не знаю как в других дистрибутивах, но у меня в Debian4.0 (дома) и в Ubuntu7.10 (на работе) команда sudo apt-get install gmt устанавливает GMT не полностью и не до конца. Данные высокого разрешения по береговым линиям и границам государств отсутствуют, и путь к GMT в системе не прописан. Поэтому, чтобы можно было нормально работать, после установки этой программы аптгетом нужно сделать следующее:

  • Прописать в системе путь к GMT. Просто вставляем в конец файла ~/.bashrc строчки:

export GMTHOME=/usr/lib/gmt/
export PATH=/usr/lib/gmt/bin/:$PATH

  • Установить данные о береговых линиях континентов высоко разрешения (взято с сache.ath.cx/Linux/map):

cd /tmp
ftp ftp.geologi.uio.no
username ftp
password guest
cd pub/gmt/3
get GMT3_high.tar.gz
get GMT3_full.tar.gz
quit
tar -xzvf GMT3_high.tar.gz
tar -xzvf GMT3_full.tar.gz
cd share
sudo mv share/*.cdf /usr/share/gmt

После этих действий GMT наконец будет установлена правильно.


Язык статистических расчетов R

25 января, 2008

Если вы желаете придать своей работе больший вес и решили для этой цели в тексте расставить в случайных местах такие слова как «корреляционный анализ», «распределение Гаусса» и «критерий Фишера», то язык для статистических расчетов R будет как раз к стати. К тому же, он умеет строить внушительные графики, при известных навыках работы с R, обескураживающие бывалых математиков. Читать далее…