Имеется:
а) совокупность однослойных растровых отображений в формате ERDAS IMG (например, модель конценрации загрязняющих веществ в приповерхностном слое атмосферы)
б) набор векторных форм (полигонов) в формате ERSI SHP (shapefile) определяющих зоны для которых нужно вычислить статистику на основе растровых отображений. Зоны определены в поле ROUTE_NAME.
Все данные имеют одну и ту же геопривязку.
Задача: вычислить статистику (MEAN, STD, MIN, MAX) для каждой из зон на основе каждого из растров. Задача не может быть выполнена вручную по причине большого числа растров и зон. Статистика должна быть сохранена в форме единой таблицы.
Решение:
1) Воспользоваться скриптом David Coley для автоматического вычисления зонной статистики для совокупности файлов (BatchZonalStatsAsTable) для ArcGIS decktop. Адаптировать его код на Python для данной задачи. Для этого изменить переменные определяющие параметры запуска скрипта выражениями вида gp.Workspace = sys.argv[3], где sys.argv[3] означает поле ввода данных 3 в интерфейсе запуска скропта в ArcGIS decktop. Скрпт запускает инструмент ArcGIS из Spatial Analyst Tool -> Zonal Statistics as Table и создаёт совокупность файлов DBF, соответствующих каждому исходному растровому образу. В таблице DBF строки соответствуют зонам в исходном SHP файле.
2) Для того чтобы объединить множество файлов DBF, находящихся в директории, в один, применяется мой скрипт на Perl. Скрипт принимает путь к директории как единсвенный параметр запуска и в стандартный вывод записывает строка за строкой значения из всех прочтённых DBF файлов, разделяя их запятыми. Запускается приблизительно так: perl dbf_merge.pl f:\working_dir > output.csv
#!/usr/local/bin/perl
dbf_merge.pl by Viacheslav Shalisko
-----------------------------------
use XBase;
my $dir = $ARGV[0];
chdir $dir || die "\nCan`t open directory $dir\n";
while (<*.dbf>) {
$file = $_;
my $database = new Xbase;
$database->open_dbf($file, undef);
my $end=$database->lastrec;
my $current_rec = 1;
while ($current_rec < $end)
{
$current_rec = $database->recno;
print "linea $current_rec de $end, ";
print "$file, ";
print &trim($database->get_field("ROUTE_NAME"));
print ", ";
print $database->get_field("COUNT");
print ", ";
print $database->get_field("AREA");
print ", ";
print $database->get_field("MIN");
print ", ";
print $database->get_field("MAX");
print ", ";
print $database->get_field("RANGE");
print ", ";
print $database->get_field("MEAN");
print ", ";
print $database->get_field("STD");
print ", ";
print $database->get_field("SUM");
print "\n";
$database->go_next;
}
$database->close_dbf;
}
sub trim {
my @out = @_;
for (@out) {
s/^\s+//;
s/\s+$//;
};
return wantarray ? @out : $out[0];
}
Читать дальше......
