Как поменять оси местами?

[Александр Мурый]

Ход процесса пока мне представляется так:

импорт SHP разреза в GRASS (в условных координатах) --> перенос координат XYZ в таблицу атрибутов с помощью v.to.db (надо сохранить полигоны, как ??) --> производим расчёты с помощью SQL (db.execute и т.п.) по формулам, что были в Excel; получаем новые координаты (географич..?) --> делаем новую область в новых координатах, экспорт туда нашего разреза (???) --> NVIZ, экспорт в VTK и т.д...

Если не сохранять полигоны, а делать всё по точкам, будет проще. Но как восстановить полигоны из точек (не v.hull же..?), я пока не представляю.

Теперь надо бы проделать это "в живую".

P.S. Я всё-таки не совсем понял: формулы в Excel -- это "приведение" к географ. коод-там, подмена осей или что? Просьба пояснить, если не трудно. А просто v.transform (афинные трансформации) не пойдёт?

[Trace]

надо сохранить полигоны, как ??

так вот как их сохранить?

Теперь надо бы проделать это "в живую".

это не сложно

P.S. Я всё-таки не совсем понял: формулы в Excel -- это "приведение" к географ. коод-там, подмена осей или что? Просьба пояснить, если не трудно. А просто v.transform (афинные трансформации) не пойдёт?

Да, это приведение к географическим координатам (проекцию писал раньше). я считаю что подмена осей - для координаты Z, так как мы ее берем из координаты Yуслов. Координаты Х и Y мы рассчитываем по формуле (линейное уравнение) из известной Xуслов.
Вот v.transform я незнаю как применить в данном случае. если вам не сложно и есть время покажите на примере.
Ексель использовал для быстрого расчета уравнения. как вы могли заметить - уравнение линии тренда. и полуавтоматического расчета файла который будет экспортироваться-импортироваться в грасс.
Возможно я Вас еще больше запутал.
Моя цель - создать объемную модель используя грасс и линукс.

[Александр Мурый]

Моя цель - создать объемную модель используя грасс и линукс.

Отличная цель. Подумаю / поверчу ещё, отпишусь, если что выйдет.

[Александр Мурый]

Кратко: я почти нашёл способ обойтись без экселя и пр. (экспериментально )
Считаю, что из этого можно даже сделать небольшой полноценный модуль для пересчёта координат через GRASS ASCII vector. Чуть доработаю и скоро выложу сюда сам скрипт.

Схема такая:
вект. карта --> v.out.ascii --> разбиваем ASCII-файл на отдел. части с помощью awk (заголовок, метаданные, данные) --> пересчитываем файл "данные" по столбцам с пом. awk --> соединяем по очереди заголовок, метаданные и новые данные в новый ASCII-файл (с пом. paste) --> импорт нового ASCII с пом. v.in.ascii --> новая вект. карта GRASS

Через SQL не пробовал, т.к. не секу.

[gamm]

сразу видно юниксера

проще (ИМХО) написать программу на R (из любого shell - это если нужно менять имена файлов. и т.д.), которая читает вектор через rgdal, правит координаты, а потом пишего его назад - и запустить R с командной строки (из того же shell). Я в качестве shell использую R

[Александр Мурый]

Я думал насчёт R. Через него, наверное, было бы реально проще.
Стыдно, но пока я не умею R вообще. А текстовый файл хотя бы прочитать можно..))

Риторический вопрос: а может, и через spgrass6 пойдёт?

Я в качестве shell использую R

Бывают случаи и потяжелее (напр., CMD.EXE)

[Александр Мурый]

Запустил R из GRASS:
library(spgrass6)
razrez <- readVECT6("razrez_clean")

Получился объект класса SpatialPolygonsDataFrame (36 полигонов, 86 границ)

Хотел получить координаты вершин для границ объекта "razrez", покопался в сети -- говорят, трудно, надо программировать.

Вопрос (видимо к gamm'у): как получить координаты вершин в R? Наверно, зависит от класса объекта, так?

[gamm]

Вопрос (видимо к gamm'у): как получить координаты вершин в R? Наверно, зависит от класса объекта, так?

во что GRASS превращается, мне неведомо, я только shape/mif читаю. Судя по документации, в такую же структуру, что и readOGR пишет

Код: Выделить всё

library(rgdal)

# --- extract polygons and run point-in-poly test
fn.mask <-"Island.shp"
lay.mask<-"Island"      #  should be the same names for some reason

# это я типа искал какие точки (subdiv.coord) попали внутрь Island

land.mask<-readOGR(fn.mask,lay.mask) # там и writeORG есть
n.poly<-length(land.mask@polygons) # там все разложено по типам, см. исходники и str(). Доков нормальных нет 
land.ind<-rep(0,nrow(subdiv.coord))
for(i in 1:n.poly) {
  n.comp<-length(land.mask@polygons[[i]]@Polygons)
  for(j in 1:n.comp) { # енто число частей в сложном объекте
    tmp.poly.xy<-land.mask@polygons[[i]]@Polygons[[j]]@coords 
    tmp.in<-(point.in.polygon(subdiv.coord[,1],subdiv.coord[,2],tmp.poly.xy[,1],tmp.poly.xy[,2]) > 0)
    if(land.mask@polygons[[i]]@Polygons[[j]]@hole) {
      land.ind[tmp.in]<-0 # point fall into the "hole" or "lake" inside land (so, not-land)
    } else {
      land.ind[tmp.in]<-1 # point fall on land
    }
  }
}

[Александр Мурый]

gamm, cпасибо, хороший повод для меня начать копаться в R..) Ещё поиграюсь с грассовскими векторами, может, пригодится дальше.

[gamm]

gamm, cпасибо, хороший повод для меня начать копаться в R..

welcome to the club

[Александр Мурый]

В продолжение темы.

Trace, я накопал по отдельности всяких приёмчиков, как считать координаты по формулам без "экселей" и "кальков". Правда, не всё в экселевском файле я понял (не дружу с матем-кой..) Так, не понял, откуда берутся коэффициенты 0.9832 и 0.18 (на которые потом умножаются значения координат).

Код: Выделить всё

# Координаты линии разреза "в плане" в системе координат, заданной в SHP-файле                                                

# x1 - коорд. крайней левой точки по оси X; y1 - коорд. крайней левой точки по оси Y                                          

# x2 - коорд. крайней правой точки по оси X; y2 - коорд. крайней правой точки по оси Y                                        

#                                                                                                                             

# Координаты разреза расчётные (x0 взять с разреза)                                                                           

# x0 - коорд. крайней левой точки по оси X; 

# y0 - коорд. крайней правой точки по оси Z; y0 = sqrt((x1-x2)^2+(y1-y2)^2)                                                    



x1=$(ogrinfo -al -so razrez_xy_m_28476.shp | grep Extent | grep -o  '[0.-9.]*' | sed 's/(//g' | sed 's/)//g' | sed -n -e '1p'\

)

x2=$(ogrinfo -al -so razrez_xy_m_28476.shp | grep Extent | grep -o  '[0.-9.]*' | sed 's/(//g' | sed 's/)//g' | sed -n -e '3p'\

)

y1=$(ogrinfo -al -so razrez_xy_m_28476.shp | grep Extent | grep -o  '[0.-9.]*' | sed 's/(//g' | sed 's/)//g' | sed -n -e '2p'\

)

y2=$(ogrinfo -al -so razrez_xy_m_28476.shp | grep Extent | grep -o  '[0.-9.]*' | sed 's/(//g' | sed 's/)//g' | sed -n -e '4p'\

)



x0=0

y0=$(echo "scale=6; sqrt(($x1-$x2)^2 + ($y1-$y2)^2)" | bc -l | xargs printf "%.6f\n")

y0 считаем с помощью консольного калькулятора bc с точностью до шестого знака после запятой (здесь). В итоге y0=146395.986513, немного отличается от полученной в экселе (146395,986772)

Вопрос: как получились коэффициенты 0.9832 и 0.18 из экселевского файла?
-----------------------------------------------------------------------------------------------------------------------------------

Поехали дальше.

v.ascii.calc.gz

(1.67 КБ) 986 скачиваний

Начальная версия скрипта (sh) для пересчёта координат по столбцам через ASCII (оно вообще надо для дела, интересно?)
Пока считает только по формулам из примера, который выше по теме. Хочу сделать ввод произвольной формулы.
С атрибутами и центроидами пока не работает, только с границами.

Код: Выделить всё

 ./v.ascii.calc -help

Description:
 Simple calculations with coordinates of vector map using ASCII files and awk

Usage:
 v.ascii.calc input=string [layer=string] output=string [--overwrite]
   [--verbose] [--quiet]

Flags:
 --o   Allow output files to overwrite existing files
 --v   Verbose module output
 --q   Quiet module output

Parameters:
   input   Input vector name
   layer   Input's layer
  output   Output vector name

Получилось странно, т.к. координаты явно разные (синие линии -- "ручной" пересчет, жёлтое -- работа скрипта):

2010-12-10-101210_1024x768_scrot.png (6.58 КБ) 14589 просмотров

Либо я с awk'ом напутал, либо что-то не учёл.. Буду смотреть ещё.

[Trace]

попробую по порядку.
1. формула которую Вы привели для расчета "y0 = sqrt((x1-x2)^2+(y1-y2)^2)", она нужна была только для расчета длины линии. Я просто не знаю как найти длину линии в Грасс. Хотя догадываюсь, что возможно через модуль v.to.db.
2.

Вопрос: как получились коэффициенты 0.9832 и 0.18 из экселевского файла?

коэффициент - a для линейного уравнения y=ax+b. Чтобы сильно "не мучаться" и быстро их рассчитать. я использовал расчет линии регрессии, если Вы посмотрите график, то уравнения приведены на нем. Калькулятор сделал мозг ленивым
для расчета этих коэффициентов мне и нужна длина линии.
для нахождения этих коэффициентов нужно решить систему линейных уравнений (для координаты Х с плана)
y1=a*x1+b
y2=a*x2+b
где х1 и х2 - координата СX2 и DX2, соответственно.
у1 - координата AX1
y2 - координата ВХ1
после небольшого преобразования мы получим b=y1-a*x1=AX1-a*СX2=AX1;
для коэффициента a выводим формулу подставив значения: a=(Y1-B)/X1
мог ошибиться. В выведении формулы

[Александр Мурый]

Trace, все написанные формулы -- это, по сути, аффинная трансформация (как описано здесь), так? Тогда мой скрипт v.ascii.calc получился велосипедом, к-й ничего нового не делает

Тогда нам бы ещё раз посмотреть в сторону модуля v.transform. Особенно параметры xshift, yshift, zshift, xscale, yscale, zscale, zrot.
Просто тема интересная и не закрыта до сих пор.

[Trace]

amuriy, мне несовсем понятно как им пользоваться и как он работает. и непонятно как ему сказать чтобы он взял координаты одной оси и подсавил с пересчетом в другую. Хотя возможно наверное это сделать используя атрибутивную таблицу. я вас правильно понимаю? но в любом случае придется расчитывать коэффициент переводной (если я не ошибаюсь).

[Александр Мурый]

Сообразил, как привязать разрез в локальных координатах к "меркатору":

Код: Выделить всё

v.transform -m in=razrez out=razrez_trans pointsfile=razrez_trans.pts

где файл трансформации

Спойлер

Код: Выделить всё

# UL NW
# UR NE
# LR SW
# LL SE
-3.17400002   1557.24499512  322917.8793 6468138.548
146413.65625 1557.24499512  466874.9068 6468138.548
-3.17400002   -77062.3125      322917.8793 6441783.519
146413.65625 -77062.3125      466874.9068 6441783.519\

Флаг "-m" выводит

Спойлер

Код: Выделить всё

CHECK MAP RESIDUALS
                Current Map                  New Map
 POINT      X coord    Y coord  |        X coord   Y coord    |      residuals

  1.         -3.17      1557.24 |    322917.88     6468138.55 |         0.00
  2.     146413.66      1557.24 |    466874.91     6468138.55 |         0.00
  3.         -3.17    -77062.31 |    322917.88     6441783.52 |         0.00
  4.     146413.66    -77062.31 |    466874.91     6441783.52 |         0.00

Спойлер

Код: Выделить всё

north=6468138.54799928
south=6441783.51900084
east=466874.90679965
west=322917.87929965
top=0.000000
bottom=0.000000

Осталось разобраться с координатами по оси Z. Как я понял, в файле трансформации коор-ты по Z не задаются, значит, надо поиграться с параметрами xscale, yscale, zscale, zrot.