Пример векторизации растра и добавления в шейп-файл полей

[Denis Rykov]

Пример векторизации растра и добавления в шейп-файл полей в пакетном режиме.

В продолжение темы Векторизация растра.

Задача
В корневой директории s:/data/MOD14A1/Baikal находятся каталоги, содержащие растровые файлы, формат имени которых AYYYY_DDD.tif (например, A2008_230.tif). Необходимо векторизовать весь этот набор данных, сохранив результат в шейп-файл, причем в шейп файл должны попасть только области растров со значениями 8 и 9. Кроме того, в шейп-файл необходимо добавить 2 поля, содержащих год и дату, которые нужно брать из имен файлов исходных растров.

Вариант решения
В ходе решения использовались библиотека GDAL 1.7.0dev, released 2008/11/26 и Python 2.5.4, установленные с помощью OSGeo4W.

Запускает командную строку и выполняем команды:

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

PATH=C:\OSGeo4W\apps\gdal-dev\bin\;%PATH%
set PYTHONPATH=C:\OSGeo4W\apps\gdal-dev\pymod

После чего запускаем файл polygonize.py, текст которого приведен ниже:

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

import csv
import os
root_path='s://data//MOD14A1//Baikal'
merged=root_path+'//'+'merged.shp'
create_merged=1
for dir in os.walk(root_path):
 for file_name in dir[2]:
  if file_name[-3:]=='tif':
    os.chdir(dir[0])
    year=file_name[1:5]
    day=str(int(file_name[6:9]))
    fname=file_name[0:-4]
    
    step1='gdal_polygonize.py'+' '+file_name+' '+ '-f "ESRI Shapefile"'+' '+fname+'_temp.shp "Active_fires" "fire_mask"'
    step2='ogr2ogr -f "ESRI Shapefile" -sql "SELECT * FROM'+' '+fname+'_temp'+' '+ 'WHERE ((fire_mask=8) or (fire_mask=9))"'+' '+fname+'_sel.shp'+' '+fname+'_temp.shp'
    step4='ogr2ogr -f CSV outputcsv'+' '+fname+'_sel.shp -lco GEOMETRY=AS_WKT'
    step5='del'+' '+fname+'_temp.shp'+' '+fname+'_temp.shx'+' '+fname+'_temp.prj'+' '+fname+'_temp.dbf'+' '+fname+'_sel.shp'+' '+fname+'_sel.shx'+' '+fname+'_sel.prj'+' '+fname+'_sel.dbf'

    os.system(step1)
    os.system(step2)
    os.system(step4)
    os.system(step5)

    os.chdir(dir[0]+'//'+'outputcsv')
    csvnamein=fname+'_sel.csv'
    csvnameout=fname+'_out.csv'
    
    f = open(csvnamein, 'rt')
    if open(csvnamein).read()!='':
     reader=csv.reader(f)
     tr=1
     f_tmp=open(csvnameout,"wt")
     for row in reader:
      if tr==1:
       f_tmp.write(row[0]+','+row[1]+',Year,Day\n');
       tr=0
      else:
       f_tmp.write('"'+row[0]+'"'+','+row[1]+','+year+','+day+'\n');
     f.close()
     f_tmp.close()
     os.system('ogr2ogr'+' '+fname+'.shp'+' '+csvnameout)
     os.system('del'+' '+csvnamein+' '+csvnameout)
     os.system('move *.* ..//')
     os.chdir('..//')
     os.system('rmdir outputcsv')
    
     if create_merged==1:
       os.system('ogr2ogr'+' '+merged+' '+fname+'.shp')
       create_merged=0
     else:
       step6='ogr2ogr -update -append'+' '+merged+' '+fname+'.shp'+' '+'-nln merged'
       os.system(step6)
    else:
     f.close()
     os.chdir('..//')
     os.system('rmdir outputcsv /q /s') 

Данный пример показывает как порой бывает удобно совместить функциональные возможности Python и консольных утилит GDAL/OGR. Если есть какие-то замечания или предложения, буду рад услышать.

P.S.: Для тех, кто хочет поэкспериментировать - один из растров во вложении.