Python+gdal/ogr vector data reading and writing __python

Common vector data formats are shapefile, Geojson, CSV, and file database gdb and spatial database PostGIS, regardless of the format of data or how to store, once open the data source, get vector layer (details refer to ogr operation vector data class structure diagram), The operation of the data is the same. The following is a detailed introduction to the reading and writing of vector data. first, open different vector files

1. Define the function that opens the data source, and iterate through all the layers, outputting their names and layers.

'
    output data source layer
    parameters: the path of the FN data source Is_write The
         mode to open the data source, 0 for read-only mode, and 1 for read-write mode
'
def print_layers ' (FN, is_write ):
    ds = Ogr. Open (FN, is_write)
    if DS is None:
        #raise OSError ("Could not open%s", fn)
        raise OSError (' could not open {} '. fo Rmat (FN)) for the
    I in range (ds. Getlayercount ()):
        Lyr = ds. Getlayer (i)
        print (' {0}:{1} '. Format (i, LYR.) GetName ()))

2, open PostGIS data source

With the above Open data source and output layer method, open the local post GIS database postgis_24_sample, the output layer is as follows:

#打开postgis数据库
    print_layers (' pg:user=postgres password=postgis dbname=postgis_24_sample ', 0)

0:tiger.county
1:tiger.state
2:tiger.place
3:tiger.cousub
4:tiger.edges
5:tiger.addrfeat
6:tiger.faces
7:tiger.zcta5
8:tiger.tract
9:tiger.tabblock
10:tiger.bg

3. Open Folder Data source (Shapefiles and CSV)

Open the Shapefiles folder, and the results of the output are as follows:

#打开shapefile文件夹
    shp_fn = R ' E:\HZZ_GXDP_Spatial_Data '
    print_layers (SHP_FN, 0)

0:communityriver
1:countyriver
2:streetriver
3:video
4:wry

second, read vector data

steps to read vector data include:

1, open the data source;

2, get the layer;

3, get the elements in the layer;

4. Get the attribute and geometry information of the element.

The following is a snippet of code that reads vector data:

Import SYS from
osgeo import ogr
import ospybook as PB

fn = R ' D:\soft\geoprocessing-with-python\china_basic _map '
ds = ogr. Open (FN, 0)
if DS is None:
    sys.exit (' could not open {0}. ') Format (FN))

' index gets the layer mode in the data source '
Lyr = ds. Getlayer (0)
print (LYR. GetName ())

i = 0 for
fea in LYR:
    pt = fea.geometry ()
    x = pt. GetX ()
    y = pt. GetY () '
    field name gets property value '
    code = FEA. GetField (' Gbcode ')
    ' object method gets the property value '
    length = fea. LENGTH
    ' object method gets the property value 2 ' '
    lpoly_ = fea[' lpoly_ ' '
    Index method gets the property value '
    Fnode_ = fea. GetField (0)
    ' Gets the property value of a particular type '
    Str_len = fea. Getfieldasstring (' length ')
    print (code, length, Fnode_, Lpoly_, Str_len, x, y)
    i + 1
    if i = +:
        Brea K

' name gets the layer in the data source '
Lyr2 = ds. Getlayer (' country ')
print (LYR2. GetName ())

Iii. Acquiring data meta data

The metadata for vector data, including the data range, geometry type, space reference, and the profile of the Layer object , and so on, the following code fragment shows how to obtain the metadata information:

Import sys from OSGeo import ogr fn = R ' D:\soft\geoprocessing-with-python\china_basic_map ' ds = ogr. Open (FN, 0) If DS is None:sys.exit (' could not open {0}. ') Format (fn)) Lyr = ds. Getlayer (0) ' Get layer range ' extent = Lyr. GetExtent () print (extent) ' (73.44696044921875, 135.08583068847656, 3.408477306365967, 53.557926177978516) ' "( Min_x, max_x, min_y, max_y) ' Get the Layer collection type, the return value is number ' ' 1: Point, 2: line, 3: face ' Geom_type = Lyr. Getgeomtype () print (Geom_type) #2 print (Geom_type = = Ogr.wkbpoint) #False print (Geom_type = = ogr.wkblinestring) #True pri NT (Geom_type = = Ogr.wkbpolygon) #False ' get geometric type through elements ' FEA = Lyr. Getfeature (0) print (Fea.geometry (). Getgeometrytype ()) #2 print (Fea.geometry (). Getgeometryname ()) #LINESTRING ' Get the space reference of the layer ' Print (LYR.      Getspatialref ()) # geogcs["gcs_wgs_1984", # datum["wgs_1984", # spheroid["wgs_84", 6378137.0,298.257223563]], #

primem["Greenwich", 0.0], # unit["Degree", 0.0174532925199433], # authority["EPSG", "4326"]] print (Lyr.schema) '' Get Layer property name and data type ' for field in Lyr.schema:print (Field.name, field. Gettypename ()) # Fnode_ Integer64 # tnode_ Integer64 # lpoly_ Integer64 # rpoly_ Integer64 # LENGTH Real # bou2_4m_ Intege R64 # bou2_4m_id Integer64 # Gbcode Integer

four, vector data writes

1, the vector data writing ideas are as follows:

(1) Open (create) a data source in read-write mode;

(2) Get the layer to add the element (or create a new layer to add the element);

(3) creating empty elements and assigning values to geometric objects and attributes;

(4) Insert the created element into the layer.

The following is a detailed code that creates a new layer using the elements in one layer:

#-*-coding:utf-8-*-import sys from osgeo import ogr ' Create a layer (created from elements in one layer) ' fn = R ' D:\soft\geoprocessing-with-pyt Hon\china_basic_map ' read-write mode to open Data source ' ds = Ogr. Open (FN, 1) If DS is None:sys.exit (' could not open {0}. ') Format (fn)) In_lyr = ds. Getlayer (' Provincial city ') # Lyr = ds. 
Getlayer (' Country ') # # from Ospybook.vectorplotter import Vectorplotter # VP = Vectorplotter (False) # Vp.plot (LYR, Fill=false) # Vp.draw () ' If there is a layer with the same name, delete the ' if ' if DS. Getlayer (' capital_city '): DS. Deletelayer (' capital_city ') Out_lyr = ds. Createlayer (' capital_city ', In_lyr. Getspatialref (), Ogr.wkbpoint) ' Creates a new layer's field ' out_lyr based on the meta-layer Field object. Createfields (In_lyr.schema) ' creates an empty element based on the layer definition ' Out_defn = out_lyr. Getlayerdefn () Out_fea = Ogr. Feature (OUT_DEFN) for in_fea in In_lyr:geom = In_fea.geometry () Out_fea. Setgeometry (Geom) for I in range (In_fea. GetFieldCount ()): value = In_fea. GetField (i) if I!= 5:out_fea. SetField (i, VAlue) Else:print (value) Out_lyr. Createfeature (OUT_FEA)

2, create a new data source

The key to creating a new data source is to use the correct driver , and each driver handles only one type of vector data. There are two ways to get the driver:(1) Get from an open dataset , which will allow you to create a new data source that is the same as the existing data source vector data format;(2) using the Getdriverbyname function in Ogr, The name of the driver that was passed to it was introduced in the OGR website. Here is the code fragment:

' Create new data source ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
import sys from
osgeo import ogr

fn = R ' D:\soft\geoprocessing-with-python\ China_basic_map '

gets the first data-driven approach: Get ' ds = OGR based on existing data sources
. Open (FN, 0)
if DS is None:
    sys.exit (' could not open {0}. ') Format (FN))

Driver = ds. Getdriver ()
print (driver. GetName ())

' Gets the second way to get a data source: Getdriverbyname '
json_driver = ogr. Getdriverbyname (' Geojson ')
print (Json_driver. GetName ())

' Create Geojson data source, the data source path should be to the specific file name '
json_fn = R ' D:\soft\geoprocessing-with-python\china_basic_ Map\json_fn.json '
json_ds = json_driver. CreateDataSource (JSON_FN)
if Json_ds is None:
    sys.exit (' could not create {0} '. Format (JSON_FN))
print ( JSON_DS)

3. New property Field

To add a field to a layer, you need a Fielddefn object that contains important information such as field name, data type, field, and precision , and the following is the code fragment for the new property field:

"' Useexceptions '
," ogr. Useexceptions ()
json_fn = R ' D:\soft\geoprocessing-with-python\china_basic_map\json_fn4.json '
json_driver = Ogr. Getdriverbyname (' Geojson ')
print (' start ')

try:
    json_ds = Json_driver. CreateDataSource (JSON_FN)
    ' new Property data field '
    LYR = Json_ds. Createlayer (' layer ')
    coord_fld = ogr. Fielddefn (' X ', ogr. oftreal)
    coord_fld. SetWidth (8)
    coord_fld. Setprecision (3)
    Lyr. CreateField (COORD_FLD)
    coord_fld. SetName (' Y ')
    lyr. CreateField (COORD_FLD)

    #新建要素
    fea = ogr. Feature (LYR. Getlayerdefn ())
    fea. SetField (' X ', 12.34)
    fea. Setfid (1)
    Lyr. Createfeature (FEA)
    #结果同步至硬盘 (saved to file)
    Json_ds. Synctodisk ()
except RuntimeError as E:
    print (e)

print (' End ')

v. Updating of existing Data

When you do vector data processing, you sometimes need to update existing data instead of creating an entirely new dataset. The ability to update and support editing operations depends on the data format. Common update data operations include changing layer definitions (property fields), factor additions, updates, and deletions . The following are described separately:

1. Change Layer definition

This includes modifying the definition of property fields, adding new property fields, and deleting property fields, as follows:

Import SYS from
osgeo import ogr
' "' Update existing Data '" '
fn = R ' D:\soft\geoprocessing-with-python\china_ Basic_map '
ds = ogr. Open (FN, 1)
if DS is None:
    sys.exit (' could not open {0}. ') Format (FN))

' Change layer definition (add, remove, modify field) '
Lyr = ds. Getlayer (0)
lyrdefn = Lyr. Getlayerdefn ()
i = Lyrdefn. Getfieldindex (' Gbcode ')
' gets the data type of the field '
Fld_type = lyrdefn. Getfielddefn (i). GetType ()
print (fld_type)
' Defines the new property field '
Fld_defn = ogr. Fielddefn (' Gbcodecode ', fld_type)
fld_defn2 = ogr. Fielddefn (' NEWFIELD4 ', ogr. Oftinteger)
' use Alter_name_flag to change the name of the property field, you must ensure that the data type of the field is consistent '
LYR. Alterfielddefn (i, Fld_defn, ogr). Alter_name_flag)
' new property field '
Lyr. CreateField (fld_defn2) '
delete field (parameter is field index value) '
lyr. DeleteField (LYR. Findfieldindex (' Newfield ', 0))
Lyr. DeleteField (Lyrdefn. Getfieldindex (' NEWFIELD2 '))

2. Factor add, UPDATE and delete

' Add, remove and update elements '
fn = R ' D:\soft\geoprocessing-with-python\china_basic_map '
ds = ogr. Open (FN, 1)
if DS is None:
    sys.exit (' could not open {0}. ') Format (FN))

#更新要素 (new property field)
Lyr = ds. Getlayer (0)
Lyr. CreateField (ogr. Fielddefn (' Newfield ', ogr. Oftinteger))
n = 1 for
fea in LYR:
    fea. SetField (' Newfield ', n*n)
    Lyr. Setfeature (FEA)
    n + + 1

#删除要素 for
fea in LYR:
    print (FEA. GetField (' Newfield '))
    if fea. GetField (' newfield ') = = 3179089:
        lyr. Deletefeature (FEA. Getfid ())

print (LYR. Getfeaturecount ())

After the data is updated, it is generally necessary to compress the database, reorganize or recalculate the data space, as follows:

Ds. ExecuteSQL (' repack ' + LYR. GetName ()) #压缩数据库
ds. ExecuteSQL (' VACUUM ') #重组功能

ds. ExecuteSQL (' Recompute EXTENT on ' + LYR. GetName ()) #重新计算图层的空间范围
print (LYR. GetExtent ())