Извлечение координат точек из растра

[ubugunov]

gornak, спасибо за совет. Пробую делать по Вашей схеме.
Растр не стал резать, решил протестировать "как есть".
Растр в точки перевел, поля добавил, вычислил геометрию, перевел в LongInteger.
Перечисленные действия выполнял в течение дня, но, полагаю, что если безотрывно сидеть ждать окончания каждой операции, то часа за два это можно сделать.
Далее начал выполнять точки в растр, часов восемь уже стоит, все закончить не может. Арка вроде не зависла, позволяет выполнять любые другие действия. Будем подождать ...

[ubugunov]

Ура, растр с координатами Lon получен!!!
Операция "Точки в растр" продлилась 13 часов.
Теперь остается создать расчетный растр искомого Hpz путем калькуляции растров SRTM и нового Lon с использованием необходимой формулы.
Предварительно прикидываю, что вся процедура на 1 растр SRTM займет минимум 2 рабочих дня. Если подгадать и поставить "точки в растр" на ночь, то, наверное, можно сократить процесс.

[gamm]

что-то вы не то делаете. Вот цепочка, которая грузит квадратик SRTM за несколько секунд, и превращает в таблицу с долготой и широтой. Посчитать, и записать результат - еще несколько секунд. Правда, это не Арка, а бесплатный R (вырезал кусок из старого скрипта, доделал - один стандартный файл SRTM 1х1 градус, 1201х1201 ячеек, обрабатывается меньше минуты)

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

library(rgdal)
image.fn<-"N39E052.hgt"

#======================================================================================
# --- GDAL stuff ----------------------------------------------------------------------
fn.region<-image.fn     # channels 1-5,7 
work.region.tif <- new("GDALReadOnlyDataset", fn.region)
p.b<-getRasterTable(work.region.tif)        # by rows, from the upper-left corner
p.bbox<-work.region.tif[,,1]@bbox           # bounding box
p.csize<-work.region.tif[,,1]@grid@cellsize # cell (pixel) size
p.b.prj<-getProjectionRef(work.region.tif)  
p.b.dim<-dim(work.region.tif)               # c(ny,nx,nband)

# --- get grid
tmp<-work.region.tif[,,1]

# --- close image
GDAL.close(work.region.tif)

print(p.bbox)
print(p.csize)
names(p.b)
# --- set proper names
names(p.b)<-c("Lon","Lat","H")

print(p.b.prj)
print(p.b.dim)


#======================================================================================
# --- Calc result=Lon+Lat+H
p.out<-data.frame(x=p.b[,1],y=p.b[,2],Res=p.b$Lon+p.b$Lat+p.b$H)

# --- make GeoTIFF image (with projection)
coordinates(p.out) = ~x+y
gridded(p.out) = TRUE
p.out.grd=as(p.out,"SpatialGridDataFrame")
p.out.grd@proj4string@projargs<-as.character(p.b.prj)
p.out.tif<-create2GDAL(p.out.grd,drivername="GTiff",type="Float32")
saveDataset(p.out.tif,"foo.tif")
GDAL.close(p.out.tif)