gdal浅尝篇

1.gdal背景知识

• GDAL主要提供了三大类数据的支持：栅格数据，矢量数据以及空间网络数据（Geographic
  Network Model）。

• GDAL不但提供了API借口方便开发人员自定义自己的功能，而且还提供了一系列实用工具（Command
  Line
  Tools）可以实现方便快速的空间数据处理。我们可以使用这些实用工具，结合Linux
  Shell脚本或者Windows批处理脚本进行大批量空间数据的批量处理。

1. gdal实现数据格式转换

• 例子：shp转换为json

• 复制需要转json的矢量文件到：并且改名为英文（例如ckx）

• 打开里面文件为：shp转json 教程.doc

  

• ogr2ogr -f “GeoJSON” ckx.json ckx.shp

• 然后在C:\Windows\System32\cmd.exe路径下打开cmd

• 运行打开修改后的这个：ogr2ogr -f “GeoJSON” ckx.json ckx.shp

• 会生成一个json文件在该路径下，然后用记事本打开，另存

• 一定改为utf-8；

• 到此使用gdal转换数据为json格式 完成。

• 其他格式数据转换

• shape 转mapinfo文件

• ogr2ogr -f “MapInfo File” -fieldTypeToString “Integer64” road Railway.shp

• shape转GeoJSON格式

• ogr2ogr -f “GeoJSON” -fieldTypeToString “Integer64” road Railway.shp

• shape转wkt格式

• ogr2ogr -lco “GEOMETRY=AS_WKT” -lco “SEPARATOR=TAB” -f CSV railway.csv -sql
  “select LUA_ID from LandUseArea” LandUseArea.shp

• Shapez转wkt格式：ogr2ogr -lco “GEOMETRY=AS_WKT” -f CSV TAZ2.csv TAZ2.shp

1. 使用arcpy代码实现shp转换为wkt

• from osgeo import ogr

• shapefile = r"shp\BOUA.shp"

• driver = ogr.GetDriverByName(“ESRI Shapefile”)

• dataSource = driver.Open(shapefile, 0)

• layer = dataSource.GetLayer()

• layerDefinition = layer.GetLayerDefn()

• #获取字段名称列表

• fields = []

• for i in range(layerDefinition.GetFieldCount()):

• fields.append(layerDefinition.GetFieldDefn(i).GetName())

• print(fields)

• #获取字段值和相应的几何

• for feature in layer:

• name = feature.GetField(“NAME”)

• pac = feature.GetField(“PAC”)

• geom = feature.GetGeometryRef()

• geomwkt = geom.ExportToWkt()

• print("%s;%s;%s\n"%(name,pac,geomwkt))

1. 使用代码将wkt转换为shp

• import osr

• import ogr

• outShapefile = r"aa4.shp"

• tabletxt = r"data.txt"

• outDriver = ogr.GetDriverByName(“ESRI Shapefile”)

• # Create the output shapefile

• outDataSource = outDriver.CreateDataSource(outShapefile)

• #设置投影

• srs = osr.SpatialReference()

• srs.ImportFromEPSG(4326)

• outLayer = outDataSource.CreateLayer(“states_extent”,
  srs,geom_type=ogr.wkbPolygon)

• # 创建3个字段并设置类型和长度

• aField = ogr.FieldDefn(“A”, ogr.OFTString)

• aField.SetWidth(5)

• outLayer.CreateField(aField)

• bField = ogr.FieldDefn(“B”, ogr.OFTString)

• bField.SetWidth(18)

• outLayer.CreateField(bField)

• cField = ogr.FieldDefn(“C”, ogr.OFTString)

• cField.SetWidth(100)

• outLayer.CreateField(cField)

• #创建几何

• featureDefn = outLayer.GetLayerDefn()

• def addPoly(a,b,c,wktpoly):

• feature = ogr.Feature(featureDefn)

• feature.SetGeometry(wktpoly)

• feature.SetField(“A”, a)

• feature.SetField(“B”, b)

• feature.SetField(“C”, c)

• outLayer.CreateFeature(feature)

• feature = None

• f = open(tabletxt,‘r’)

• #文件可能比较大，逐行读取

• line = f.readline()

• while line:

• a,b,c,wktstr = line.split(";")

• wktpoly = ogr.CreateGeometryFromWkt(wktstr)

• addPoly(a,b,c, wktpoly)

• line = f.readline()

• f.close()

• outDataSource = None

1. 采用arcgis中字段计算器来计算wkt参数

打开shp文件的字段——添加文本类型字段——字段长度为4000——计算窗口示意如下:

同理还可计算：!SHAPE.JSON!