from osgeo import gdal
from osgeo import ogr
# 导入ogr,注意导入的名称
from osgeo import ogr
# 载入GeoJSON驱动
driver = 0.GetDriverByName('GeoJSON')
print(driver)
<osgeo.ogr.Driver; proxy of <Swig Object of type 'OGRDriverShadow *' at 0x1149c3c90> >
# 不分大小写,这样写也行
driver = ogr.GetDriverByName('geojson')
print(driver)
<osgeo.ogr.Driver; proxy of <Swig Object of type 'OGRDriverShadow *' at 0x1149c3db0> >
driver2 = ogr.GetDriverByName('Esri shapefile')
print(driver2)
<osgeo.ogr.Driver; proxy of <Swig Object of type 'OGRDriverShadow *' at 0x1149c3cc0> >
# 从shapefile中读取点要素,以城市人口数据为例
import sys
from osgeo import ogr
# 打开数据源,打开图层
fn = 'E:/公众号/python/2/osgeopy-data/global/ne_50m_populated_places.shp'
# ogr.Open函数,打开数据源,注意Open的O是大写
# 根据文件类型,Open函数自动载入驱动
# 第一个参数是shapefile的文件名
# 第二个参数update =0:只读,=1:读写
# 如果文件无法打开,ds is None,注意错误处理
# 如果文件打开正确,返回数据源ds(data source)
ds = ogr.Open(fn, 0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
# 一个数据源包括多个图层,用GetLayer(n)函数返回索引为n的图层,该参数默认值为0
lyr = ds.GetLayer(0)
i = 0
# 循环读取每个要素
for feat in lyr:
# 读取x和y方向的坐标
pt = feat.geometry()
x = pt.GetX()
y = pt.GetY()
# 读取属性值
name = feat.GetField('NAME')
pop = feat.GetField('POP_MAX')
print(name, pop, x, y)
i += 1
if i == 10:
break
del ds
Bombo 75000 32.533299524864844 0.5832991056146284
Fort Portal 42670 30.27500161597942 0.671004121125236
Potenza 69060 15.798996495640267 40.642002130098206
Campobasso 50762 14.655996558921856 41.56299911864397
Aosta 34062 7.315002595706176 45.7370010670723
Mariehamn 10682 19.949004471869102 60.09699618489543
Ramallah 24599 35.206209378189556 31.90294475142406
Vatican City 832 12.453386544971766 41.903282179960115
Poitier 85960 0.3332765285345545 46.58329225573658
Clermont-Ferrand 233050 3.080008095928406 45.779982115759424
# 打开shp文件
data_dir = 'E:/公众号/python/2/osgeopy-data/'
fn = os.path.join(data_dir, 'global', 'ne_50m_poplated_places.shp')
# 根据文件类型,Open函数会自动载入驱动
ds = ogr.Open(fn, 0)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
lyr = ds.GetLayer(0)
# Get the total number of features and the last one.
# 载入要素的总数
num_features = lyr.GetFeatureCount()
print(num_features)
# 获取最后一个要素, 按索引值访问要素
last_feature = lyr.GetFeature(num_features - 1)
print(last_feature.NAME)
1249
Hong Kong
# 注意:会记得最后访问的要素对象(称为当前要素)下次读取会从当前要素之后开始读
# 使用lyr.ResetReading()重置当前要素
# 如果循环两遍,却没有重置,则第二遍不会打印任何东西
fn = os.path.join(data_dir, 'Washington', 'large_cities.geojson')
ds = ogr.Open(fn, 0)
lyr = ds.GetLayer(0)
print('First loop')
for feat in lyr:
print(feat.GetField('Name'), feat.GetField('Population'))
print('Second loop')
for feat in lyr:
pt = feat.geometry()
print(feat.GetField('Name'), pt.GetX(), pt.GetY())
First loop
Bellingham 67171
Everett 91488
Shoreline 53025
Seattle 563374
Bellevue 109569
Renton 50052
Kent 79524
Federal Way 83259
Tacoma 193556
Spokane 195629
Yakima 71845
Vancouver 143560
Kennewick 54693
Second loop
# 重置之后,正常执行下一个循环
ds = ogr.Open(fn, 0)
lyr = ds.GetLayer(0)
print('First loop')
for feat in lyr:
print(feat.GetField('Name'), feat.GetField('Population'))
print('Second loop')
lyr.ResetReading() # This is the important line.
for feat in lyr:
pt = feat.geometry()
print(feat.GetField('Name'), pt.GetX(), pt.GetY())
First loop
Bellingham 67171
Everett 91488
Shoreline 53025
Seattle 563374
Bellevue 109569
Renton 50052
Kent 79524
Federal Way 83259
Tacoma 193556
Spokane 195629
Yakima 71845
Vancouver 143560
Kennewick 54693
Second loop
Bellingham -122.48822021484375 48.759552001953125
Everett -122.20207214355469 47.97898864746094
Shoreline -122.34152221679688 47.75564956665039
Seattle -122.33206939697266 47.6062126159668
Bellevue -122.20067596435547 47.610382080078125
Renton -122.2170639038086 47.482879638671875
Kent -122.2348403930664 47.38093185424805
Federal Way -122.31262969970703 47.32231903076172
Tacoma -122.44429779052734 47.25288009643555
Spokane -117.4260482788086 47.658782958984375
Yakima -120.50589752197266 46.60206604003906
Vancouver -122.66148376464844 45.638729095458984
Kennewick -119.1372299194336 46.211246490478516
# Open the large_cities data source.
fn = os.path.join(data_dir, 'Washington', 'large_cities.geojson')
ds = ogr.Open(fn)
if ds is None:
sys.exit('Could not open {0}.'.format(fn))
# 读取空间范围
lyr = ds.GetLayer(0)
extent = lyr.GetExtent()
print(extent)
print('Upper left corner: {}, {}'.format(extent[0], extent[3]))
print('Lower right corner: {}, {}'.format(extent[1], extent[2]))
(-122.66148376464844, -117.4260482788086, 45.638729095458984, 48.759552001953125)
Upper left corner: -122.66148376464844, 48.759552001953125
Lower right corner: -117.4260482788086, 45.638729095458984
# 读取几何要素类型
print(lyr.GetGeomType())
print(lyr.GetGeomType() == ogr.wkbPoint)
print(lyr.GetGeomType() == ogr.wkbPolygon)
1
True
False
# 几何要素类型,易读模式
feat = lyr.GetFeature(0)
print(feat.geometry().GetGeometryName())
POINT
# 读取空间参考
print(lyr.GetSpatialRef())
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4269"]]
# 读取字段名称和类型
for field in lyr.schema:
print(field.name, field.GetTypeName())
CITIESX020 Integer
FEATURE String
NAME String
POP_RANGE String
POPULATION Integer
FIPS55 String
COUNTY String
FIPS String
STATE String
STATE_FIPS String
DISPLAY Integer
ds.CreateLayer(name,[srs],[geom_type],[options])
# 将一个shp文件中的要素写入另一个文件
import sys
from osgeo import ogr
# 打开数据源,注意update参数设置为1,表示可读可写
# 这一次的数据源是文件夹,不是文件,在该文件夹下创建新的shp文件
ds = ogr.Open(data_dir+'/global', 1)
if ds is None:
sys.exit('Could not open folder.')
# 将已有的shp文件视作图层,读取
in_lyr = ds.GetLayer('ne_50m_populated_places')
# 新建图层,如果同名图层已经存在,先把它删掉
if ds.GetLayer('capital_cities'):
ds.DeleteLayer('capital_cities')
# 以后使用图层时,建议引用图层名称而不是索引,写入新图层后,索引会变化
# 新建图层,设定图层名称,空间参考,类型
out_lyr = ds.CreateLayer('capital_cities',
in_lyr.GetSpatialRef(),
ogr.wkbPoint)
# 新建图层的属性字段定义列表,将其设定为和源图层一模一样
out_lyr.CreateFields(in_lyr.schema)
# 新建空的要素
out_defn = out_lyr.GetLayerDefn()
out_feat = ogr.Feature(out_defn)
for in_feat in in_lyr:
# 如果FEATURECLA属性为‘Admin-0 capital’ 则进行以下操作
if in_feat.GetField('FEATURECLA') == 'Admin-0 capital':
# Copy geometry and attributes
# 复制几何信息
geom = in_feat.geometry()
out_feat.SetGeometry(geom)
# 复制属性信息
for i in range(in_feat.GetFieldCount()):
value = in_feat.GetField(i)
out_feat.SetField(i, value)
# 插入要素CreateFeature()
out_lyr.CreateFeature(out_feat)
# 插入要素之后,该要素的信息就会写入图层对象。out_feat对象可以重复使用,不会影响已经写入的信息。
# 最后需要保证写入的信息同步到文件中
# 强制关闭数据源:
del ds
# 或者使用同步函数,不关闭数据源
# 一定要及时的同步,否则会造成写入
# 创建一个新的数据源,而不是写入到已有的数据源
# 从已有数据源提取驱动程序
ds = ogr.Open(os.path.join(data_dir, 'global', 'ne_50m_admin_0_countries.shp'))
driver = ds.GetDriver()
print(driver.name)
ESRI Shapefile
# 载入指定名称的驱动程序
json_driver = ogr.GetDriverByName('GeoJSON')
print(json_driver.name)
# 创建GeoJSON数据源
# 注意这个文件夹必须要首先存在,ogr不会自动创建文件夹
json_fn = os.path.join(data_dir, 'output', 'example.geoson')
json_ds = json_driver.CreateDataSource(json_fn)
if json_ds is None:
sys.exit('Could not create {0}.'.format(json_fn))
GeoJSON
# 如果数据源已经存在,先删掉它
if os.path.exists(json_fn):
json_driver.DeleteDataSource(json_fn)
# 创建SQLite数据库
driver = ogr.GetDriverByName('SQLite')
ds = driver.CreateDataSource(os.path.join(data_dir, 'out-put', 'earth.sqlite'),
['SPATIALITE=yes'])
# 使用ogr的抛出异常功能
ogr.UseExceptions(
fn = os.path.join(data_dir, 'output', 'africa.geojson')
driver = ogr.GetDriverByName('GeoJSON')
print('Doing some preliminary analysis...')
try:
# 如果文件已经存在,会抛出异常
ds = driver.CreateDataSource(fn)
lyr = ds.CreateLayer('layer')
except RuntimeError as e:
# 在这里处理异常
print(e)
print('Doing some more analysis...')
Doing some preliminary analysis...
Doing some more analysis...
ogr.FieldDefn()
coord_fld = ogr.FieldDefn('X_short', ogr.OFTReal)
coord_fld.SetXXXX()
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
lyr.CreateField(fld)
out_lyr.CreateField(coord_fld)
coord_fld.SetName('Y_short')
out_lyr.CreateField(coord_fld)
# 打开输入的shp文件,创建输出的shp文件,和以前类似
in_fn = os.path.join(data_dir, 'Washington', 'large_cities.shp')
in_ds = ogr.Open(in_fn, 0)
if in_ds is None:
sys.exit('Could not open {0}.'.format(in_fn))
in_lyr = in_ds.GetLayer(0)
driver = in_ds.GetDriver()
out_fn = os.path.join(data_dir, 'output', 'precision_test.shp')
if os.path.exists(out_fn):
driver.DeleteDataSource(out_fn)
out_ds = driver.CreateDataSource(out_fn)
if out_ds is None:
sys.exit('Could not create {0}.'.format(out_fn))
out_lyr = out_ds.CreateLayer('precision_test',
in_lyr.GetSpatialRef(),
ogr.wkbPoint)
# 使用FieldDefn创建字段定义
name_fld = ogr.FieldDefn('Name', ogr.OFTString)
# 使用SetXXXX系列函数改变字段的设置
name_fld.SetWidth(6)
# 使用CreateField创建字段
out_lyr.CreateField(name_fld)
# 用默认的精度创建两个字段,这两个字段只是名字不一样,用SetName改字段定义的名字
coord_fld = ogr.FieldDefn('X_default', ogr.OFTReal)
out_lyr.CreateField(coord_fld)
coord_fld.SetName('Y_default')
out_lyr.CreateField(coord_fld)
# 创建两个字段,修改精度(其实是降低了精度)
coord_fld = ogr.FieldDefn('X_short', ogr.OFTReal)
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
out_lyr.CreateField(coord_fld)
coord_fld.SetName('Y_short')
out_lyr.CreateField(coord_fld)
# 将数据复制进去。用QGIS或者ArcGIS打开文件,看看里面的内容
out_feat = ogr.Feature(out_lyr.GetLayerDefn())
for in_feat in in_lyr:
pt = in_feat.geometry()
name = in_feat.GetField('NAME')
out_feat.SetGeometry(in_feat.geometry())
out_feat.SetField('Name', name)
out_feat.SetField('X_default', pt.GetX())
out_feat.SetField('Y_default', pt.GetY())
out_feat.SetField('X_short', pt.GetX())
out_feat.SetField('Y_short', pt.GetY())
out_lyr.CreateFeature(out_feat)
lyr.AlterFieldDefn(i, fld_defn, ogr.FLAG)
original_fn = os.path.join(data_dir, 'Washington', 'large_cities.shp')
new_fn = os.path.join(data_dir, 'output', 'large_cities2.shp')
# 将原始文件复制一份
pb.copy_datasource(original_fn, new_fn)
ds = ogr.Open(new_fn, 1)
if ds is None:
sys.exit('Could not open {0}.'.format(new_fn))
lyr = ds.GetLayer(0)
# 查看一下属性
print('Original attributes')
pb.print_attributes(lyr, geom=True)
Original attributes
FID Geometry CITIESX020 FEATURE NAME POP_RANGE POP_2000 FIPS55 COUNTY FIPS STATE STATE_FIPS DISPLAY
0 POINT (-122.488, 48.760) 404 County Seat Bellingham 50,000 - 99,999 67171 05280 Whatcom County 53073 WA 53 0
1 POINT (-122.202, 47.979) 538 County Seat Everett 50,000 - 99,999 91488 22640 Snohomish County 53061 WA 53 0
2 POINT (-122.342, 47.756) 613 Populated Place Shoreline 50,000 - 99,999 53025 63960 King County 53033 WA 53 0
3 POINT (-122.332, 47.606) 685 County Seat Seattle 500,000 - 999,999 563374 63000 King County 53033 WA 53 0
4 POINT (-122.201, 47.610) 697 Populated Place Bellevue 100,000 - 249,999 109569 05210 King County 53033 WA 53 0
5 POINT (-122.217, 47.483) 776 Populated Place Renton 50,000 - 99,999 50052 57745 King County 53033 WA 53 0
6 POINT (-122.235, 47.381) 838 Populated Place Kent 50,000 - 99,999 79524 35415 King County 53033 WA 53 0
7 POINT (-122.313, 47.322) 870 Populated Place Federal Way 50,000 - 99,999 83259 23515 King County 53033 WA 53 0
8 POINT (-122.444, 47.253) 896 County Seat Tacoma 100,000 - 249,999 193556 70000 Pierce County 53053 WA 53 0
9 POINT (-117.426, 47.659) 1184 County Seat Spokane 100,000 - 249,999 195629 67000 Spokane County 53063 WA 53 0
10 POINT (-120.506, 46.602) 1594 County Seat Yakima 50,000 - 99,999 71845 80010 Yakima County 53077 WA 53 0
11 POINT (-122.661, 45.639) 2035 County Seat Vancouver 100,000 - 249,999 143560 74060 Clark County 53011 WA 53 0
12 POINT (-119.137, 46.211) 2103 Populated Place Kennewick 50,000 - 99,999 54693 35275 Benton County 53005 WA 53 0
13 of 13 features
# 新建一个字段定义,使用AlterFieldDefn函数修改已存在的图层中的字段
# 将字段名从Name改为City_Name
i = lyr.GetLayerDefn().GetFieldIndex('Name')
fld_defn = ogr.FieldDefn('City_Name', ogr.OFTString)
lyr.AlterFieldDefn(i, fld_defn, ogr.ALTER_NAME_FLAG)
# 查看一下属性
print('Modified attributes')
pb.print_attributes(lyr, geom=True)
Modified attributes
FID Geometry CITIESX020 FEATURE City_Name POP_RANGE POP_2000 FIPS55 COUNTY FIPS STATE STATE_FIPS DISPLAY
0 POINT (-122.488, 48.760) 404 County Seat Bellingham 50,000 - 99,999 67171 05280 Whatcom County 53073 WA 53 0
1 POINT (-122.202, 47.979) 538 County Seat Everett 50,000 - 99,999 91488 22640 Snohomish County 53061 WA 53 0
2 POINT (-122.342, 47.756) 613 Populated Place Shoreline 50,000 - 99,999 53025 63960 King County 53033 WA 53 0
3 POINT (-122.332, 47.606) 685 County Seat Seattle 500,000 - 999,999 563374 63000 King County 53033 WA 53 0
4 POINT (-122.201, 47.610) 697 Populated Place Bellevue 100,000 - 249,999 109569 05210 King County 53033 WA 53 0
5 POINT (-122.217, 47.483) 776 Populated Place Renton 50,000 - 99,999 50052 57745 King County 53033 WA 53 0
6 POINT (-122.235, 47.381) 838 Populated Place Kent 50,000 - 99,999 79524 35415 King County 53033 WA 53 0
7 POINT (-122.313, 47.322) 870 Populated Place Federal Way 50,000 - 99,999 83259 23515 King County 53033 WA 53 0
8 POINT (-122.444, 47.253) 896 County Seat Tacoma 100,000 - 249,999 193556 70000 Pierce County 53053 WA 53 0
9 POINT (-117.426, 47.659) 1184 County Seat Spokane 100,000 - 249,999 195629 67000 Spokane County 53063 WA 53 0
10 POINT (-120.506, 46.602) 1594 County Seat Yakima 50,000 - 99,999 71845 80010 Yakima County 53077 WA 53 0
11 POINT (-122.661, 45.639) 2035 County Seat Vancouver 100,000 - 249,999 143560 74060 Clark County 53011 WA 53 0
12 POINT (-119.137, 46.211) 2103 Populated Place Kennewick 50,000 - 99,999 54693 35275 Benton County 53005 WA 53 0
13 of 13 features
# 同时改变字段的多个属性(例如同时更改名称和精度)
lyr_defn = lyr.GetLayerDefn()
i = lyr_defn.GetFieldIndex('POP_2000')
width = lyr_defn.GetFieldDefn(i).GetWidth()
fld_defn = ogr.FieldDefn('POP_2000x', ogr.OFTReal)
fld_defn.SetWidth(width)
fld_defn.SetPrecision(4)
flag = ogr.ALTER_NAME_FLAG + ogr.ALTER_WIDTH_PRECISION_FLAG
lyr.AlterFieldDefn(i, fld_defn, flag)
# 查看一下属性
print('Modified attributes')
pb.print_attributes(lyr, geom=True)
Modified attributes
FID Geometry CITIESX020 FEATURE City_Name POP_RANGE POP_2000x FIPS55 COUNTY FIPS STATE STATE_FIPS DISPLAY
0 POINT (-122.488, 48.760) 404 County Seat Bellingham 50,000 - 99,999 67171 05280 Whatcom County 53073 WA 53 0
1 POINT (-122.202, 47.979) 538 County Seat Everett 50,000 - 99,999 91488 22640 Snohomish County 53061 WA 53 0
2 POINT (-122.342, 47.756) 613 Populated Place Shoreline 50,000 - 99,999 53025 63960 King County 53033 WA 53 0
3 POINT (-122.332, 47.606) 685 County Seat Seattle 500,000 - 999,999 563374 63000 King County 53033 WA 53 0
4 POINT (-122.201, 47.610) 697 Populated Place Bellevue 100,000 - 249,999 109569 05210 King County 53033 WA 53 0
5 POINT (-122.217, 47.483) 776 Populated Place Renton 50,000 - 99,999 50052 57745 King County 53033 WA 53 0
6 POINT (-122.235, 47.381) 838 Populated Place Kent 50,000 - 99,999 79524 35415 King County 53033 WA 53 0
7 POINT (-122.313, 47.322) 870 Populated Place Federal Way 50,000 - 99,999 83259 23515 King County 53033 WA 53 0
8 POINT (-122.444, 47.253) 896 County Seat Tacoma 100,000 - 249,999 193556 70000 Pierce County 53053 WA 53 0
9 POINT (-117.426, 47.659) 1184 County Seat Spokane 100,000 - 249,999 195629 67000 Spokane County 53063 WA 53 0
10 POINT (-120.506, 46.602) 1594 County Seat Yakima 50,000 - 99,999 71845 80010 Yakima County 53077 WA 53 0
11 POINT (-122.661, 45.639) 2035 County Seat Vancouver 100,000 - 249,999 143560 74060 Clark County 53011 WA 53 0
12 POINT (-119.137, 46.211) 2103 Populated Place Kennewick 50,000 - 99,999 54693 35275 Benton County 53005 WA 53 0
13 of 13 features
# 使用copy方法拷贝字段定义,可以避免一些奇怪的错误
import copy
lyr_defn = lyr.GetLayerDefn()
i = lyr_defn.GetFieldIndex('POP_2000x')
fld_defn = copy.copy(lyr_defn.GetFieldDefn(i))
fld_defn.SetName('POP_2000xx')
fld_defn.SetPrecision(4)
flag = ogr.ALTER_NAME_FLAG + ogr.ALTER_WIDTH_PRECISION_FLAG
lyr.AlterFieldDefn(i, fld_defn, flag)
print('\nNew field names and empty ID field')
pb.print_attributes(lyr, geom=False)
New field names and empty ID field
FID CITIESX020 FEATURE City_Name POP_RANGE POP_2000xx FIPS55 COUNTY FIPS STATE STATE_FIPS DISPLAY
0 404 County Seat Bellingham 50,000 - 99,999 67171 05280 Whatcom County 53073 WA 53 0
1 538 County Seat Everett 50,000 - 99,999 91488 22640 Snohomish County 53061 WA 53 0
2 613 Populated Place Shoreline 50,000 - 99,999 53025 63960 King County 53033 WA 53 0
3 685 County Seat Seattle 500,000 - 999,999 563374 63000 King County 53033 WA 53 0
4 697 Populated Place Bellevue 100,000 - 249,999 109569 05210 King County 53033 WA 53 0
5 776 Populated Place Renton 50,000 - 99,999 50052 57745 King County 53033 WA 53 0
6 838 Populated Place Kent 50,000 - 99,999 79524 35415 King County 53033 WA 53 0
7 870 Populated Place Federal Way 50,000 - 99,999 83259 23515 King County 53033 WA 53 0
8 896 County Seat Tacoma 100,000 - 249,999 193556 70000 Pierce County 53053 WA 53 0
9 1184 County Seat Spokane 100,000 - 249,999 195629 67000 Spokane County 53063 WA 53 0
10 1594 County Seat Yakima 50,000 - 99,999 71845 80010 Yakima County 53077 WA 53 0
11 2035 County Seat Vancouver 100,000 - 249,999 143560 74060 Clark County 53011 WA 53 0
12 2103 Populated Place Kennewick 50,000 - 99,999 54693 35275 Benton County 53005 WA 53 0
13 of 13 features
# 增加要素
lyr.ResetReading()
lyr.CreateField(ogr.FieldDefn('ID', ogr.OFTInteger))
n = 1
for feat in lyr:
feat.SetField('ID', n)
lyr.SetFeature(feat)
n += 1
print('\nID has been added and precision has taken effect')
pb.print_attributes(lyr, geom=False)
ID has been added and precision has taken effect
FID CITIESX020 FEATURE City_Name POP_RANGE POP_2000xx FIPS55 COUNTY FIPS STATE STATE_FIPS DISPLAY ID
0 404 County Seat Bellingham 50,000 - 99,999 67171 05280 Whatcom County 53073 WA 53 0 1
1 538 County Seat Everett 50,000 - 99,999 91488 22640 Snohomish County 53061 WA 53 0 2
2 613 Populated Place Shoreline 50,000 - 99,999 53025 63960 King County 53033 WA 53 0 3
3 685 County Seat Seattle 500,000 - 999,999 563374 63000 King County 53033 WA 53 0 4
4 697 Populated Place Bellevue 100,000 - 249,999 109569 05210 King County 53033 WA 53 0 5
5 776 Populated Place Renton 50,000 - 99,999 50052 57745 King County 53033 WA 53 0 6
6 838 Populated Place Kent 50,000 - 99,999 79524 35415 King County 53033 WA 53 0 7
7 870 Populated Place Federal Way 50,000 - 99,999 83259 23515 King County 53033 WA 53 0 8
8 896 County Seat Tacoma 100,000 - 249,999 193556 70000 Pierce County 53053 WA 53 0 9
9 1184 County Seat Spokane 100,000 - 249,999 195629 67000 Spokane County 53063 WA 53 0 10
10 1594 County Seat Yakima 50,000 - 99,999 71845 80010 Yakima County 53077 WA 53 0 11
11 2035 County Seat Vancouver 100,000 - 249,999 143560 74060 Clark County 53011 WA 53 0 12
12 2103 Populated Place Kennewick 50,000 - 99,999 54693 35275 Benton County 53005 WA 53 0 13
13 of 13 features
# 删除要素,例如删除名称为西雅图的城市
lyr.ResetReading()
for feat in lyr:
if feat.GetField('City_Name') == 'Seattle':
lyr.DeleteFeature(feat.GetFID())
print('\nSeattle deleted')
pb.print_attributes(lyr, geom=False)
# 删除之后需要整理/重组数据库空间,并重新计算空间范围
ds.ExecuteSQL('REPACK ' + lyr.GetName())
ds.ExecuteSQL('RECOMPUTE EXTENT ON ' + lyr.GetName())
print('\nDatabase packed')
pb.print_attributes(lyr, geom=False)
Seattle deleted
FID CITIESX020 FEATURE City_Name POP_RANGE POP_2000xx FIPS55 COUNTY FIPS STATE STATE_FIPS DISPLAY ID
0 404 County Seat Bellingham 50,000 - 99,999 67171 05280 Whatcom County 53073 WA 53 0 1
1 538 County Seat Everett 50,000 - 99,999 91488 22640 Snohomish County 53061 WA 53 0 2
2 613 Populated Place Shoreline 50,000 - 99,999 53025 63960 King County 53033 WA 53 0 3
3 697 Populated Place Bellevue 100,000 - 249,999 109569 05210 King County 53033 WA 53 0 5
4 776 Populated Place Renton 50,000 - 99,999 50052 57745 King County 53033 WA 53 0 6
5 838 Populated Place Kent 50,000 - 99,999 79524 35415 King County 53033 WA 53 0 7
6 870 Populated Place Federal Way 50,000 - 99,999 83259 23515 King County 53033 WA 53 0 8
7 896 County Seat Tacoma 100,000 - 249,999 193556 70000 Pierce County 53053 WA 53 0 9
8 1184 County Seat Spokane 100,000 - 249,999 195629 67000 Spokane County 53063 WA 53 0 10
9 1594 County Seat Yakima 50,000 - 99,999 71845 80010 Yakima County 53077 WA 53 0 11
10 2035 County Seat Vancouver 100,000 - 249,999 143560 74060 Clark County 53011 WA 53 0 12
11 2103 Populated Place Kennewick 50,000 - 99,999 54693 35275 Benton County 53005 WA 53 0 13
12 of 12 features
Database packed
FID CITIESX020 FEATURE City_Name POP_RANGE POP_2000xx FIPS55 COUNTY FIPS STATE STATE_FIPS DISPLAY ID
0 404 County Seat Bellingham 50,000 - 99,999 67171 05280 Whatcom County 53073 WA 53 0 1
1 538 County Seat Everett 50,000 - 99,999 91488 22640 Snohomish County 53061 WA 53 0 2
2 613 Populated Place Shoreline 50,000 - 99,999 53025 63960 King County 53033 WA 53 0 3
3 697 Populated Place Bellevue 100,000 - 249,999 109569 05210 King County 53033 WA 53 0 5
4 776 Populated Place Renton 50,000 - 99,999 50052 57745 King County 53033 WA 53 0 6
5 838 Populated Place Kent 50,000 - 99,999 79524 35415 King County 53033 WA 53 0 7
6 870 Populated Place Federal Way 50,000 - 99,999 83259 23515 King County 53033 WA 53 0 8
7 896 County Seat Tacoma 100,000 - 249,999 193556 70000 Pierce County 53053 WA 53 0 9
8 1184 County Seat Spokane 100,000 - 249,999 195629 67000 Spokane County 53063 WA 53 0 10
9 1594 County Seat Yakima 50,000 - 99,999 71845 80010 Yakima County 53077 WA 53 0 11
10 2035 County Seat Vancouver 100,000 - 249,999 143560 74060 Clark County 53011 WA 53 0 12
11 2103 Populated Place Kennewick 50,000 - 99,999 54693 35275 Benton County 53005 WA 53 0 13
12 of 12 features
ogrinfo myfile.xxx
ogrinfo -so myfile.xxx
ogrinfo -fid 1 -q -geom=No myfile.xxx
# 过滤数据举例,先画个图
ds = ogr.Open(os.path.join(data_dir, 'global'))
lyr = ds.GetLayer('ne_50m_admin_0_countries')
vp.plot(lyr, fill=False)
pb.print_attributes(lyr, 4, ['name'], geom=False)
FID name
0 Aruba
1 Afghanistan
2 Angola
3 Anguilla
4 of 241 features
# 按要素属性字段过滤,选择亚洲国家
# lyr.SetAttributeFilter('continent = "Asia"') # wrong!
lyr.SetAttributeFilter("continent = 'Asia'") # right!
lyr.GetFeatureCount()
# Draw the Asian countries in yellow and print out a few features.
vp.plot(lyr, 'y')
pb.print_attributes(lyr, 4, ['name'], geom=False)
FID name
1 Afghanistan
7 United Arab Emirates
9 Armenia
17 Azerbaijan
4 of 53 features
# 其他要素都还在,你可以按照FID也就是要素编号选择
print(lyr.GetFeature(2).GetField('name'))
Angola
# 新建一个过滤器,这次选南美洲
# lyr.SetAttributeFilter('continent = "South America"') # wrong!
lyr.SetAttributeFilter("continent = 'South America'") # right!
vp.plot(lyr, 'b')
lyr.SetAttributeFilter(None)
lyr.GetFeatureCount()
241
# 空间过滤
# 打开数据源,载入数据
ds = ogr.Open(os.path.join(data_dir, 'global'))
country_lyr = ds.GetLayer('ne_50m_admin_0_countries')
vp.plot(country_lyr, fill=False)
# 绘制世界城市
city_lyr = ds.GetLayer('ne_50m_populated_places')
city_lyr.GetFeatureCount()
vp.plot(city_lyr, 'y.')
# 创建过滤器,过滤德国
country_lyr.SetAttributeFilter("name = 'Germany'")
feat = country_lyr.GetNextFeature()
germany = feat.geometry().Clone()
# 从德国的多边形创建过滤器
city_lyr.SetSpatialFilter(germany)
print(city_lyr.GetFeatureCount())
vp.plot(country_lyr, fill=False)
vp.plot(city_lyr, 'bo')
5
# 再创建一个过滤器,只选择百万人口以上的城市
city_lyr.SetAttributeFilter('pop_min > 1000000')
print(city_lyr.GetFeatureCount())
vp.plot(country_lyr, fill=False)
vp.plot(city_lyr, 'rs')
3
# 去掉过滤器
city_lyr.SetSpatialFilter(None)
city_lyr.SetAttributeFilter(None)
print(city_lyr.GetFeatureCount())
country_lyr.SetSpatialFilter(None)
country_lyr.SetAttributeFilter(None)
print(country_lyr.GetFeatureCount())
vp.plot(country_lyr, fill=False)
vp.plot(city_lyr, 'm^', markersize=8)
1249
241
# 按经纬度范围过滤
# 先清除所有过滤器
vp.clear() # 清空绘图
country_lyr.SetAttributeFilter(None)
vp.plot(country_lyr, fill=False)
# 设置经纬度范围作为过滤器
country_lyr.SetSpatialFilterRect(110, -50, 160, 10)
vp.plot(country_lyr, 'y')
ExecuteSQL(statement, [spatialFilter], [dialect])
# 创建点要素“火堆”,使用ogr.Geometry函数,设定要素类型为“点”
firepit = ogr.Geometry(ogr.wkbPoint)
# 设定点的坐标
# 向点对象中添加一个点的坐标
# 点对象只有一个点坐标
# 修改点对象也同样使用AddPoint,原有坐标被覆盖
firepit.AddPoint(59.5, 11.5)
# 打印出来看看
print(firepit)
# 可视化
vp = VectorPlotter(True)
vp.plot(firepit, 'bo')
POINT (59.5 11.5 0)
# 修改点的坐标也是用AddPoint函数,原有坐标被覆盖
firepit.AddPoint(59.5, 13)
vp.plot(firepit, 'rs')
print(firepit)
POINT (59.5 13.0 0)
# 还可以使用SetPoint函数编辑点的坐标
firepit.SetPoint(0, 59.5, 33)
print(firepit)
POINT (59.5 33.0 0)
# 构造2.5D点对象,加入高程信息
firepit = ogr.Geometry(ogr.wkbPoint25D)
firepit.AddPoint(59.5, 11.5, 2)
print(firepit)
POINT (59.5 11.5 2)
Geometry type | 2D constant | 2.5D constant |
---|---|---|
Point | wkbPoint | wkbPoint25D |
Multipoint | wkbMultiPoint | wkbMultiPoint25D |
Line | wkbLineString | wkbLineString25D |
Multiline | wkbMultiLineString | wkbMultiLineString25D |
Polygon ring | wkbLinearRing | n/a |
Polygon | wkbPolygon | wkbPolygon25D |
Multipolygon | wkbMultiPolygon | wkbMultiPolygon25D |
Geometry collection | wkbGeometryCollection | wkbGeometryCollection25D |
# 创建多点要素“水龙头”,由多个单点要素合并而成
faucets = ogr.Geometry(ogr.wkbMultiPoint)
faucet = ogr.Geometry(ogr.wkbPoint)
faucet.AddPoint(67.5, 16)
faucets.AddGeometry(faucet)
faucet.AddPoint(73, 31)
faucets.AddGeometry(faucet)
faucet.AddPoint(91, 24.5)
faucets.AddGeometry(faucet)
0
# 可视化
vp.clear()
vp.plot(faucets, 'bo')
vp.zoom(-5)
print(faucets)
MULTIPOINT (67.5 16.0 0,73 31 0,91.0 24.5 0)
# 更改其中一个点的坐标
faucets.GetGeometryRef(1).AddPoint(25, 32)
vp.plot(faucets, 'k^', 'tmp')
print(faucets)
MULTIPOINT (67.5 16.0 0,25 32 0,91.0 24.5 0)
# 批量更改坐标
faucets2 = faucets.Clone()
for i in range(faucets2.GetGeometryCount()):
pt = faucets2.GetGeometryRef(i)
pt.AddPoint(pt.GetX() + 2, pt.GetY())
# 比较更改前后的坐标
vp.plot(faucets, 'bs')
vp.plot(faucets2, 'rs')
vp.zoom(-5)
# 创建线要素“人行道”,注意添加点的顺序不能乱
sidewalk = ogr.Geometry(ogr.wkbLineString)
sidewalk.AddPoint(54, 37)
sidewalk.AddPoint(62, 35.5)
sidewalk.AddPoint(70.5, 38)
sidewalk.AddPoint(74.5, 41.5)
# 可视化
vp = VectorPlotter(True)
vp.plot(sidewalk, 'b-')
print(sidewalk)
LINESTRING (54 37 0,62.0 35.5 0,70.5 38.0 0,74.5 41.5 0)
# 更改点的坐标(更改的是3号点)
sidewalk.SetPoint(3, 74.5, 31.5)
vp.plot(sidewalk, 'k--', 'tmp')
print(sidewalk)
LINESTRING (54 37 0,62.0 35.5 0,70.5 38.0 0,74.5 31.5 0)
# 批量修改线要素中每个点的坐标
# 使用SetPoint指定点的索引,逐个修改
sidewalk2 = sidewalk.Clone()
for i in range(sidewalk.GetPointCount()):
sidewalk.SetPoint(i, sidewalk.GetX(i), sidewalk.GetY(i)+ 1)
vp.plot(sidewalk, 'b--')
vp.plot(sidewalk2, 'r--')
print(sidewalk)
LINESTRING (54 38 0,62.0 36.5 0,70.5 39.0 0,74.5 32.5 0)
# 清点要素中点的个数
print(sidewalk.GetPointCount()) # vertices
# 清点要素中sub-geometries(子几何对象)的个数
print(sidewalk.GetGeometryCount()) # sub-geometries
4
0
# 将线对象的顶点导出为列表
print(sidewalk.GetPoints())
[(54.0, 38.0, 0.0), (62.0, 36.5, 0.0), (70.5, 39.0, 0.0), (74.5, 32.5, 0.0)]
# 在第二和第三个顶点之间插入一个点
# 注意插入的语法是[2:2]
vertices = sidewalk.GetPoints()
vertices[2:3] = [(66.5, 35)]
print(vertices)
[(54.0, 38.0, 0.0), (62.0, 36.5, 0.0), (66.5, 35), (74.5, 32.5, 0.0)]
# 从上面的列表生成新的线对象
new_sidewalk = ogr.Geometry(ogr.wkbLineString)
for vertex in vertices:
new_sidewalk.AddPoint(*vertex)
vp.plot(new_sidewalk, 'g:')
# 从geojson文件中读取线对象
# 注意!!!原geojson文件语法有错误,倒数第三行最后多了个一逗号
ds = ogr.Open(os.path.join(data_dir, 'misc', 'line-example.geojson'))
lyr = ds.GetLayer()
feature = lyr.GetFeature(0)
line = feature.geometry().Clone()
vp.clear()
vp.plot(line, 'b-')
# 为线对象添加新的顶点
# 使用GetPoints得到顶点列表,用列表的方式添加
# 用*符号把传递的元组对象“炸开”
# 从后往前添加,因为加入新顶点后,后边的顶点索引会自动更新
# 如果你从前往后加,容易搞乱了到底是往哪个点后面添加新的点
vertices = line.GetPoints()
vertices[26:26] = [(87, 57)]
vertices[19:19] = [(95, 38), (97, 43), (101, 42)]
vertices[11:11] = [(121, 18)]
vertices[5:5] = [(67, 32), (74, 30)]
new_line = ogr.Geometry(ogr.wkbLineString)
for vertex in vertices:
new_line.AddPoint(*vertex)
vp.plot(line, 'b-')
vp.plot(new_line, 'b--')
# 创建多个Line对象,然后合并为MultiLine
path1 = ogr.Geometry(ogr.wkbLineString)
path1.AddPoint(61.5, 29)
path1.AddPoint(63, 20)
path1.AddPoint(62.5, 16)
path1.AddPoint(60, 13)
path2 = ogr.Geometry(ogr.wkbLineString)
path2.AddPoint(60.5, 12)
path2.AddPoint(68.5, 13.5)
path3 = ogr.Geometry(ogr.wkbLineString)
path3.AddPoint(69.5, 33)
path3.AddPoint(80, 33)
path3.AddPoint(86.5, 22.5)
paths = ogr.Geometry(ogr.wkbMultiLineString)
paths.AddGeometry(path1)
paths.AddGeometry(path2)
paths.AddGeometry(path3)
0
# 可视化
vp.clear()
vp.plot(paths, 'b-')
print(paths)
MULTILINESTRING ((61.5 29.0 0,63 20 0,62.5 16.0 0,60 13 0),(60.5 12.0 0,68.5 13.5 0),(69.5 33.0 0,80 33 0,86.5 22.5 0))
# 多线对象的编辑
# 获取多线对象内部的线对象,然后逐一编辑
# 例如,获取多线内部的0号对象(0号线),然后更改索引为1的点的坐标(从0开始数)
paths2 = paths.Clone()
paths2.GetGeometryRef(0).SetPoint(1, 73, 28)
vp.plot(paths, 'b-')
vp.plot(paths2, 'k--', 'tmp')
print(paths)
print(paths2)
MULTILINESTRING ((61.5 29.0 0,63 20 0,62.5 16.0 0,60 13 0),(60.5 12.0 0,68.5 13.5 0),(69.5 33.0 0,80 33 0,86.5 22.5 0))
MULTILINESTRING ((61.5 29.0 0,73 28 0,62.5 16.0 0,60 13 0),(60.5 12.0 0,68.5 13.5 0),(69.5 33.0 0,80 33 0,86.5 22.5 0))
# 批量编辑点的坐标
paths3 = paths.Clone()
for i in range(paths3.GetGeometryCount()):
path = paths3.GetGeometryRef(i)
for j in range(path.GetPointCount()):
path.SetPoint(j, path.GetX(j) + 2, path.GetY(j) - 3)
vp.plot(paths, 'b-')
vp.plot(paths3, 'r--')
# 首先新建一个环,添加顶点
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(58, 38.5)
ring.AddPoint(53, 6)
ring.AddPoint(99.5, 19)
ring.AddPoint(73, 42)
# 新建一个多边形
yard = ogr.Geometry(ogr.wkbPolygon)
# 将环添加到多边形
yard.AddGeometry(ring)
# 将环闭合
yard.CloseRings()
# 可视化一下
vp = VectorPlotter(False)
vp.plot(yard, fill=False, edgecolor='blue')
print(yard)
POLYGON ((58.0 38.5 0,53 6 0,99.5 19.0 0,73 42 0,58.0 38.5 0))
# 编辑多边形对象
# 要先取出环再编辑
# 环不用放回去,更新了环之后争夺个多边形也更新了
yard2 = yard.Clone()
ring = yard2.GetGeometryRef(0)
for i in range(ring.GetPointCount()):
ring.SetPoint(i, ring.GetX(i) - 5, ring.GetY(i))
vp.plot(yard, fill=False, edgecolor='blue')
vp.plot(yard2, fill=False, ec='red', linestyle='dashed')
# 向多边形添加新的顶点,先取出环,再编辑
yard3 = yard.Clone()
ring = yard3.GetGeometryRef(0)
vertices = ring.GetPoints()
vertices[2:3] = ((90, 16), (90, 27))
for i in range(len(vertices)):
ring.SetPoint(i, *vertices[i])
vp.plot(yard, fill=False, edgecolor='blue')
vp.plot(yard3, fill=False, ec='black', ls='dotted', linewidth=3)
# 新建复合多边形对象
# 甜甜圈状的嵌套多边形
# 创建两个环,先外环再内环,然后闭合
# 第一个环
box1 = ogr.Geometry(ogr.wkbLinearRing)
box1.AddPoint(87.5, 25.5)
box1.AddPoint(89, 25.5)
box1.AddPoint(89, 24)
box1.AddPoint(87.5, 24)
garden1 = ogr.Geometry(ogr.wkbPolygon)
garden1.AddGeometry(box1)
# 第二个环
box2 = ogr.Geometry(ogr.wkbLinearRing)
box2.AddPoint(89, 23)
box2.AddPoint(92, 23)
box2.AddPoint(92,22)
box2.AddPoint(89,22)
garden2 = ogr.Geometry(ogr.wkbPolygon)
garden2.AddGeometry(box2)
# 闭合
gardens = ogr.Geometry(ogr.wkbMultiPolygon)
gardens.AddGeometry(garden1)
gardens.AddGeometry(garden2)
gardens.CloseRings()
vp.clear()
vp.plot(gardens, fill=False, ec='blue')
vp.zoom(-1)
print(gardens)
MULTIPOLYGON (((87.5 25.5 0,89.0 25.5 0,89 24 0,87.5 24.0 0,87.5 25.5 0)),((89 23 0,92 23 0,92 22 0,89 22 0,89 23 0)))
# 编辑复合多边形,要取出多边形,再从中取出每个环,进行编辑
gardens2 = gardens.Clone()
for i in range(gardens2.GetGeometryCount()):
ring = gardens2.GetGeometryRef(i).GetGeometryRef(0)
for j in range(ring.GetPointCount()):
ring.SetPoint(j, ring.GetX(j) + 1, ring.GetY(j) + 0.5)
vp.plot(gardens, fill=False, ec='blue')
vp.plot(gardens2, fill=False, ec='red', ls='dashed')
vp.zoom(-1)
# 甜甜圈状多边形的编辑
# 用for循环取出每一个环,逐个编辑
# 不要假定你的多边形不是甜甜圈,这会在后期带来bug
lot = ogr.Geometry(ogr.wkbLinearRing)
lot.AddPoint(58, 38.5)
lot.AddPoint(53, 6)
lot.AddPoint(99.5, 19)
lot.AddPoint(73, 42)
house = ogr.Geometry(ogr.wkbLinearRing)
house.AddPoint(67.5, 29)
house.AddPoint(69, 25.5)
house.AddPoint(64, 23)
house.AddPoint(69, 15)
house.AddPoint(82.5, 22)
house.AddPoint(76, 31.5)
yard = ogr.Geometry(ogr.wkbPolygon)
yard.AddGeometry(lot)
yard.AddGeometry(house)
yard.CloseRings()
vp.clear()
vp.plot(yard, 'yellow')
print(yard)
yard2 = yard.Clone()
for i in range(yard2.GetGeometryCount()):
ring = yard2.GetGeometryRef(i)
for j in range(ring.GetPointCount()):
ring.SetPoint(j, ring.GetX(j) - 5, ring.GetY(j))
print(yard2)
vp.plot(yard2, fill=False, hatch='x', color='blue')
POLYGON ((58.0 38.5 0,53 6 0,99.5 19.0 0,73 42 0,58.0 38.5 0),(67.5 29.0 0,69.0 25.5 0,64 23 0,69 15 0,82.5 22.0 0,76.0 31.5 0,67.5 29.0 0))
POLYGON ((53.0 38.5 0,48 6 0,94.5 19.0 0,68 42 0,53.0 38.5 0),(62.5 29.0 0,64.0 25.5 0,59 23 0,64 15 0,77.5 22.0 0,71.0 31.5 0,62.5 29.0 0))
* 相交:Intersects
* 接触:Touches
* 穿过:Crosses
* 内部:Within
* 包含:Contains
* 重叠:Overlaps
* 不相交:Disjoint
# 载入新奥尔良周边的湿地边界(已经预先裁剪好)
vp = VectorPlotter(False)
water_ds = ogr.Open(os.path.join(data_dir, 'US', 'wtrbdyp010.shp'))
water_lyr = water_ds.GetLayer(0)
water_lyr.SetAttributeFilter('WaterbdyID = 1011327')
marsh_feat = water_lyr.GetNextFeature()
marsh_geom = marsh_feat.geometry().Clone()
vp.plot(marsh_geom, 'b')
# 载入新奥尔良市行政区范围边界
nola_ds = ogr.Open(os.path.join(data_dir, 'Louisiana', 'NOLA.shp'))
nola_lyr = nola_ds.GetLayer(0)
nola_feat = nola_lyr.GetNextFeature()
nola_geom = nola_feat.geometry().Clone()
vp.plot(nola_geom, fill=False, ec='red', ls='dashed', lw=3)
# 使用相交分析函数
# 获得新奥尔良城市内的湿地范围
# 注意面向对象的语法
intersection = marsh_geom.Intersection(nola_geom)
vp.plot(marsh_geom, 'b')
vp.plot(nola_geom, fill=False, ec='red', ls='dashed', lw=3)
vp.plot(intersection, 'yellow', hatch='x')
# 统计新奥尔良市行政区范围内有多少是湿地
# 使用GetArea方法,通过空间过滤器实现
water_lyr.SetAttributeFilter("Feature != 'Lake'")
water_lyr.SetSpatialFilter(nola_geom)
wetlands_area = 0
for feat in water_lyr:
intersect = feat.geometry().Intersection(nola_geom)
wetlands_area += intersect.GetArea()
print(intersect.GetArea())
pcnt = wetlands_area / nola_geom.GetArea()
print('{:.1%} of New Orleans is wetland'.format(pcnt))
0.0006515340173894044
0.0014098862841060988
0.022241193829142178
6.680527991115466e-07
28.7% of New Orleans is wetland
# 统计新奥尔良市行政区范围内有多少是湿地
# 在内存中建立临时图层,使用SQL语句实现
water_lyr.SetSpatialFilter(None)
water_lyr.SetAttributeFilter("Feature != 'Lake'")
memory_driver = ogr.GetDriverByName('Memory')
temp_ds = memory_driver.CreateDataSource('temp')
temp_lyr = temp_ds.CreateLayer('temp')
nola_lyr.Intersection(water_lyr, temp_lyr)
sql = 'SELECT SUM(OGR_GEOM_AREA) AS area FROM temp'
lyr = temp_ds.ExecuteSQL(sql)
pcnt = lyr.GetFeature(0).GetField('area') / nola_geom.GetArea()
print('{:.1%} of New Orleans is wetland'.format(pcnt))
28.7% of New Orleans is wetland
# 美国有多少城市在火山周边10英里(16000米)内?
# 一种有缺陷的方法
# 每个火山周围画出16000米的缓冲区,然后做叠加分析
from osgeo import ogr
# 载入数据
shp_ds = ogr.Open(os.path.join(data_dir,'US'))
volcano_lyr = shp_ds.GetLayer('us_volcanos_albers')
cities_lyr = shp_ds.GetLayer('cities_albers')
# 在内存中建立一个新的图层
# Create a temporary layer in memory to store the buffers.
memory_driver = ogr.GetDriverByName('memory')
memory_ds = memory_driver.CreateDataSource('temp')
buff_lyr = memory_ds.CreateLayer('buffer')
buff_feat = ogr.Feature(buff_lyr.GetLayerDefn())
# 对所有的火山循环一遍,每个火山周围建立半径为16000米的缓冲区
for volcano_feat in volcano_lyr:
buff_geom = volcano_feat.geometry().Buffer(16000)
tmp = buff_feat.SetGeometry(buff_geom)
tmp = buff_lyr.CreateFeature(buff_feat)
# 建好缓冲区后,进行叠加分析,城市(点图层)和火山周边区域(多边形图层)相交的点,进行计数
result_lyr = memory_ds.CreateLayer('result')
buff_lyr.Intersection(cities_lyr, result_lyr)
print('Cities: {}'.format(result_lyr.GetFeatureCount()))
Cities: 83
# 在上一个例子中,有的城市离多个火山的距离不超过10英里,被算了多次!
# 美国有多少城市在火山周边10英里(16000米)内?
# 一种更好的方法
# 每个火山周围画出16000米的缓冲区
# 然后合并为一个复合多边形对象,再做叠加分析
import ogr
shp_ds = ogr.Open(os.path.join(data_dir,'US'))
volcano_lyr = shp_ds.GetLayer('us_volcanos_albers')
cities_lyr = shp_ds.GetLayer('cities_albers')
multipoly = ogr.Geometry(ogr.wkbMultiPolygon)
for volcano_feat in volcano_lyr:
buff_geom = volcano_feat.geometry().Buffer(16000)
multipoly.AddGeometry(buff_geom)
cities_lyr.SetSpatialFilter(multipoly.UnionCascaded())
print('Cities: {}'.format(cities_lyr.GetFeatureCount()))
Cities: 78
34° 01’ 43.72995’’ -117° 12’ 54.61539’’
34° 01’ 43.77884’’ -117° 12’ 57.75961’’
# 载入数据,读取SRS
ds = ogr.Open(os.path.join(data_dir, 'US', 'states_48.shp'))
srs = ds.GetLayer().GetSpatialRef()
# 打印SRS,WKT格式
print(srs)
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4269"]]
# 用不同的格式打印
# PROJ.4
print(srs.ExportToProj4())
+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
# 创建空间参考对象
from osgeo import osr
utm_sr = osr.SpatialReference()
utm_sr.ImportFromEPSG(26912)
print(utm_sr) # WKT
print(utm_sr.ExportToProj4()) # PROJ.4
PROJCS["NAD83 / UTM zone 12N",
# 打印空间参考的相关属性
# Get the projection name.
print(utm_sr.GetAttrValue('PROJCS'))
# Get the authority.
print(utm_sr.GetAttrValue('AUTHORITY'))
print(utm_sr.GetAttrValue('AUTHORITY', 1))
# Get the datum code.
print(utm_sr.GetAuthorityCode('DATUM'))
NAD83 / UTM zone 12N
EPSG
26912
6269
# Create a UTM SRS from an EPSG code.
sr = osr.SpatialReference()
sr.ImportFromEPSG(26912)
print(sr.GetAttrValue('PROJCS'))
NAD83 / UTM zone 12N
# Create a UTM SRS from a PROJ.4 string.
sr = osr.SpatialReference()
sr.ImportFromProj4('''+proj=utm +zone=12 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs ''')
print(sr.GetAttrValue('PROJCS'))
UTM Zone 12, Northern Hemisphere
# Create a unprojected SRS from a WKT string.
wkt = '''GEOGCS["GCS_North_American_1983",
DATUM["North_American_Datum_1983",
SPHEROID["GRS_1980",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]]'''
sr = osr.SpatialReference(wkt)
print(sr)
GEOGCS["GCS_North_American_1983",
DATUM["North_American_Datum_1983",
SPHEROID["GRS_1980",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]]
# Create an Albers SRS using parameters.
sr = osr.SpatialReference()
sr.SetProjCS('USGS Albers')
sr.SetWellKnownGeogCS('NAD83')
sr.SetACEA(29.5, 45.5, 23, -96, 0, 0)
sr.Fixup()
sr.Validate()
print(sr)
PROJCS["USGS Albers",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4269"]],
PROJECTION["Albers_Conic_Equal_Area"],
PARAMETER["standard_parallel_1",29.5],
PARAMETER["standard_parallel_2",45.5],
PARAMETER["latitude_of_center",23],
PARAMETER["longitude_of_center",-96],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
# 绘制全球陆地
world = pb.get_shp_geom(os.path.join(data_dir, 'global', 'ne_110m_land_1p.shp'))
vp = VectorPlotter(True)
vp.plot(world)
# 埃菲尔铁塔的坐标转换,设定经纬度坐标
tower = ogr.Geometry(wkt='POINT (2.294694 48.858093)')
tower.AssignSpatialReference(osr.SpatialReference(osr.SRS_WKT_WGS84))
# 按投影坐标系,Web墨卡托投影绘制
# 注意:全球和埃菲尔铁搭这两个几何对象本身都被改变了
gdal.SetConfigOption('OGR_ENABLE_PARTIAL_REPROJECTION','TRUE')
web_mercator_sr = osr.SpatialReference()
web_mercator_sr.ImportFromEPSG(3857)
world = pb.get_shp_geom(os.path.join(data_dir, 'global','ne_110m_land_1p.shp'))
world.TransformTo(web_mercator_sr)
tower.TransformTo(web_mercator_sr)
vp.clear()
vp.plot(world)
print(tower)
# Web墨卡托其实不是一个正式的投影坐标系,
# 是Google码农任性开脑洞的产物,所以它的代码也比较的有意思…
# 测绘专家和码农们,你们先撕,我们吃瓜群众看看就好
#(参考 https://blog.csdn.net/kikitamoon/article/details/46124935 )
POINT (255444.16760638 6250816.95768025)
# 地理坐标转换到投影坐标
import pyproj
utm_proj = pyproj.Proj('+proj=utm +zone=31 +ellps=WGS84')
x, y = utm_proj(2.294694, 48.858093)
print(x, y)
448265.9146797116 5411920.651461346
# 投影坐标转换到地理坐标
# 注意将inverse选项设置为True
x1, y1 = utm_proj(x, y, inverse=True)
print(x1, y1)
2.294694 48.85809299999998
# 不同坐标系之间的转换:WGS84转换为NAD27
wgs84 = pyproj.Proj('+proj=utm +zone=18 +datum=WGS84')
nad27 = pyproj.Proj('+proj=utm +zone=18 +datum=NAD27')
x, y = pyproj.transform(wgs84, nad27, 580744.32, 4504695.26)
print(x, y)
580712.7534284436 4504476.041497709
今天的分享的关于利用python处理矢量数据就到这里啦!内容有些多,希望对大家有帮助。