Составление простого растра по очень большой БД

[biogis]

Добрый вечер!
Имеется БД местонахождений биообъектов из iNaturalist, содержащая порядка 60 миллионов записей в формате CSV, включающих и координатную информацию о местоположениях биообъектов – http://www.inaturalist.org/observations ... s-dwca.zip - На основании этих данных надо составить глобальный растр в общегеографической проекции с пространственным разрешением 30 секунд, формат geotif. В клетки растра надо ввести суммарное число местонахождений биообъектов, приходящихся на соответствующую клетку растра.
Вроде бы простая задача, но объем записей огромен. Пока не понятно: в каком софте эту БД открывать и обрабатывать с целью получения искомого растра.

[gamm]

1) у вас такой раст мало куда влезет, а если и влезет, то работать с ним тяжело будет, надо резать на тайлы (квадратики), SRTM (3-сек), например, нарезается градус на градус, т.е. 1200х1200, вам можно 10х10 градусов, например. Поэтому вопрос номер один - что вы с ним делать собираетесь?

2) очевидно, что таблицу нужно импортировать в любую современную базу, арифметикой посчитать номера строк и столбцов, а потом запросом посчитать для них количества. Postgress точно потянет, скорее всего и MySQL тоже, а то и Access. Потом арифметикой считаем номера тайлов, и в цикле выдаем xyz файл для каждого. Можно просто программу написать на С++, и сделать то же самое за один проход, раскидав xyz по файлам.

3) Читаем xyz из п.2 в C++/R/Python, и создаем растровые тайлы.

[biogis]

Спасибо!
1) итоговый глобальный растр, на самом деле, будет очень небольшой. Разрешение 30 секунд - то есть примерно 1-километровый. Проблема только в присвоении ячейкам значений количества приходящихся на них местонахождений (observations).
Растр нужен для нормализации и межрегионального сравнения данных о распространении отдельных видов с учетом активности корреспондентов БД в каждом регионе. Этот суммарный растр плотностей наблюдений как раз и должен служить показателем региональной активности корреспондентов. Пример использования: в Мичигане по изучаемому виду 100 наблюдений, в провинции Гуанчжоу - 10. Вопрос: где встречаемость вида выше? Допустим, растр показывает, что плотность регистрации прочих объектов на единицу площади в Мичигане в 50 раз выше, чем в Гуанчжоу - то есть активность корреспондентов в 50 раз выше. Тогда можно примерно считать, что встречаемость объекта в Гуаньчжоу в 5 раз выше, чем в Мичигане: 10/(100/50)=5. - Примерно так.
2) Access попробовал использовать - он такую БД не открывает.
3) при 30-секундном разрешении тайлы не потребуются - всё уложится в один растр.

[trir]

PostGIS
https://postgis.net/docs/using_raster_d ... ng_Rasters

[gamm]

сделал на скорую руку на R, но он память кушает очень. Пришлось сохранять промежуточный результат, и перезапускать, даже 64Гб не хватает. Грузилось примерно полчаса, остальное минуты. В одну клетку что-то многовато попало, а так вроде ничего. На разные детали типа вида не смотрел, все свалено в одну кучу.

промежуточный
tiff

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

fn="I:/temp2/observations.csv"
beg.time=Sys.time()
p=read.csv(fn,header=TRUE,sep=",",quote = "\"",dec = ".",fill=TRUE,comment.char="",as.is=TRUE)
end.time=Sys.time(); print(end.time-beg.time)
p1=p[,c("decimalLatitude","decimalLongitude")]

fn2="I:/temp2/observations.RData"
save(p,file=fn2)
rm(p)
gc()

p2=as.matrix(p1)
m1=cbind(p2,p2,p2[,1])
colnames(m1)[3:5]=c("Col","Row","One")
step=1/(60*2)
m1[,"Col"]=round((m1[,"decimalLongitude"]+180)/step)
m1[,"Row"]=round((m1[,"decimalLatitude"] + 90)/step)
m1[,"One"]=1
sel=!is.na(m1[,"decimalLongitude"]) & !is.na(m1[,"decimalLatitude"])
m2=m1[sel,]

beg.time=Sys.time()
m3=aggregate(m2[,"One"],by=list(Col=m2[,"Col"],Row=m2[,"Row"]),length)
end.time=Sys.time(); print(end.time-beg.time)

fn="I:/temp2/observations_counts.RData"
save(m3,file=fn)

=======================================================================
library(raster)

fn="I:/temp2/observations_counts.RData"
load(file=fn)

r1=raster(nrows=21600,ncols=43200,xmn=-180,xmx=180,ymn=-90,ymx=90)
r1[]=0

col=pmax(1,m3[,"Col"])
row=pmax(1,m3[,"Row"])
cnt=m3[,"x"]
pos=(row-1)*ncol(r1)+col
beg.time=Sys.time()
r1[pos]=cnt
end.time=Sys.time(); print(end.time-beg.time)
dataType(r1)="INT4S"

fn.out="I:/temp2/observations.tif"
beg.time=Sys.time()
tmp=writeRaster(r1,filename=fn.out,format="GTiff",overwrite=TRUE)
end.time=Sys.time(); print(end.time-beg.time)

А в целом уважаемый trir как всегда прав, вся сила в SQL

[biogis]

Спасибо!) Я, к сожалению, с этим софтом не работаю.
На пробнике север с югом для точек поменялись местами – получилось зеркально) То что все виды (корреспонденты) суммируются - так и надо. У меня при экспорте растра в Idrisi в сумме получилось 38636978.063869 (real) наблюдений – меньше, чем в БД. Да, заметил - максимальное число наблюдений на клетку составило 37831 – очень много, но, возможно, что так и есть в каких-то крупных городах - надо смотреть.
И да, - надо осваивать софт под большие данные. Или дружить)

[gamm]

насчет количества - около 100,000 строк было без координат (меньше, но точно не помню). Насчет зеркального отражения - вполне может быть, я заполнял как матрицу, у нее нумерация зеркальная.

[trir]

надо осваивать софт под большие данные

60 миллионов записей - это не большие данные, обычная СУБД спокойно справляется

[biogis]

Насчет зеркального отражения - вполне может быть, я заполнял как матрицу, у нее нумерация зеркальная

Транспонировал пробный растр. Попробую использовать. О результатах сообщу, если интересно. Еще раз спасибо!)

[tikhpetr]

gamm, Убили Вы меня. Промозглой февральской ночью финкой в переулке. Вся бездна незнания вскрылась моего. Что можно почитать? Не сочтите за труд!

[gamm]

попробовал немного поправить код - более правильно считать номера ячеек, и рпзвернул растр по вертикали. Посмотрите, если что-то еще надо, пишите. Интересно, по старой памяти, я же тоже типа эколог.
промежуточный
растр
Посмотрел картинку, вроде теперь Америка на месте. Вообще, похоже что много данных есть только по Северной Америке, Европе, Австралии и Японии, остальное просто белые пятна.

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

fn="I:/temp2/observations.csv"
beg.time=Sys.time()
p=read.csv(fn,header=TRUE,sep=",",quote = "\"",dec = ".",fill=TRUE,comment.char="",as.is=TRUE)
end.time=Sys.time(); print(end.time-beg.time)
p1=p[,c("decimalLatitude","decimalLongitude")]

fn2="I:/temp2/observations.RData"
save(p,file=fn2)
rm(p)
gc()

p2=as.matrix(p1)
m1=cbind(p2,p2,p2[,1])
colnames(m1)[3:5]=c("Col","Row","One")
step=1/(60*2)
m1[,"Col"]=floor((m1[,"decimalLongitude"]+180)/step)+1
m1[,"Row"]=floor((m1[,"decimalLatitude"] + 90)/step)+1
m1[,"One"]=1
sel=!is.na(m1[,"decimalLongitude"]) & !is.na(m1[,"decimalLatitude"])
m2=m1[sel,]

beg.time=Sys.time()
m3=aggregate(m2[,"One"],by=list(Col=m2[,"Col"],Row=m2[,"Row"]),length)[quote=tikhpetr post_id=183684 time=1677116414 user_id=17276]Что можно почитать?[/quote]
end.time=Sys.time(); print(end.time-beg.time)

fn="I:/temp2/observations_counts.RData"
save(m3,file=fn)

=======================================================================
library(raster)

fn="I:/temp2/observations_counts.RData"
load(file=fn)

r1=raster(nrows=21600,ncols=43200,xmn=-180,xmx=180,ymn=-90,ymx=90)
r1[]=0

col=pmax(1,m3[,"Col"])
row=pmax(1,m3[,"Row"])
cnt=m3[,"x"]
pos=(nrow(r1)-row+1)*ncol(r1)+col
beg.time=Sys.time()
r1[pos]=cnt
end.time=Sys.time(); print(end.time-beg.time)
dataType(r1)="INT4S"

fn.out="I:/temp2/observations.tif"
beg.time=Sys.time()
tmp=writeRaster(r1,filename=fn.out,format="GTiff",overwrite=TRUE)
end.time=Sys.time(); print(end.time-beg.time)

Что можно почитать?

я читал в основном чужие исхоники, и писал свои. Но сейчас это уже никому не нужно, за редким исключением, всем нужна готовая кнопка (и в 95% случаев она действительно где-то есть). Потому такие знания уже не нужны, так что не горюйте

[biogis]

Посмотрите

Спасибо! - качаю подкорректированный растр.

похоже что много данных есть только по Северной Америке, Европе, Австралии и Японии

Так и есть. Ещё Южная Корея. В Африке в ЮАР побольше, чем в других странах континента. - Все эти различия следует учитывать для нормализации данных и составления выборки в ЭГ анализе. Для чего и нужен растр плотностей наблюдений. Проверю растр на своих объектах и сообщу. Как с Вами связаться?

такие знания уже не нужны

Нужны) новые продвинутые кнопки строятся на новых знаниях и подходах, а познание бесконечно

[gamm]

Как с Вами связаться?

Хотел написать, но на Вашем сайте почты не нашел. Пишите сюда glmmivanov@yandex.ru

P.S. Знакомые эксерсизы, sampling efforts нужен оказался

[tikhpetr]

gamm, Вам Бог завез с математикой. Мне со слухом. Игорь и Эдуард, извините за флейм. По состоянию на февраль 2023 года я могу опознать около 30-40 почерков. Только это нахрен никому не нать.

[tikhpetr]

З.Ы. Сам пока работаю без подчерка.