zoukankan      html  css  js  c++  java
  • 投影变换

    #!/usr/bin/env python
    # -*- coding: utf-8 -*-

    import os

    try:
    import ogr
    import osr

    except ImportError:
    from osgeo import ogr, osr


    shp_driver = ogr.GetDriverByName('ESRI Shapefile')

    # input SpatialReference
    input_srs = osr.SpatialReference()
    input_srs.ImportFromEPSG(4326)

    # output SpatialReference
    output_srs = osr.SpatialReference()
    output_srs.ImportFromEPSG(3857)

    # create the CoordinateTransformation
    coord_trans = osr.CoordinateTransformation(input_srs, output_srs)

    # get the input layer
    input_shp = shp_driver.Open(r'../geodata/UTM_Zone_Boundaries.shp')
    in_shp_layer = input_shp.GetLayer()

    # create the output layer
    output_shp_file = r'../geodata/UTM_Zone_Boundaries_3857.shp'
    # check if output file exists if yes delete it
    if os.path.exists(output_shp_file):
    shp_driver.DeleteDataSource(output_shp_file)

    # create a new Shapefile object
    output_shp_dataset = shp_driver.CreateDataSource(output_shp_file)

    # create a new layer in output Shapefile and define its geometry type
    output_shp_layer = output_shp_dataset.CreateLayer("UTM_Zone_Boundaries_3857", geom_type=ogr.wkbMultiPolygon)

    # add fields to the new output Shapefile
    # get list of attribute fields
    # create new fields for output
    in_layer_def = in_shp_layer.GetLayerDefn()
    for i in range(0, in_layer_def.GetFieldCount()):
    field_def = in_layer_def.GetFieldDefn(i)
    output_shp_layer.CreateField(field_def)

    # get the output layer's feature definition
    output_layer_def = output_shp_layer.GetLayerDefn()

    # loop through the input features
    in_feature = in_shp_layer.GetNextFeature()
    while in_feature:
    # get the input geometry
    geom = in_feature.GetGeometryRef()
    # reproject the geometry
    geom.Transform(coord_trans)
    # create a new feature
    output_feature = ogr.Feature(output_layer_def)
    # set the geometry and attribute
    output_feature.SetGeometry(geom)
    for i in range(0, output_layer_def.GetFieldCount()):
    output_feature.SetField(output_layer_def.GetFieldDefn(i).GetNameRef(), in_feature.GetField(i))
    # add the feature to the shapefile
    output_shp_layer.CreateFeature(output_feature)
    # destroy the features and get the next input feature
    output_feature.Destroy()
    in_feature.Destroy()
    in_feature = in_shp_layer.GetNextFeature()

    # close the shapefiles
    input_shp.Destroy()
    output_shp_dataset.Destroy()

    spatialRef = osr.SpatialReference()
    spatialRef.ImportFromEPSG(3857)

    spatialRef.MorphToESRI()
    prj_file = open('UTM_Zone_Boundaries.prj', 'w')
    prj_file.write(spatialRef.ExportToWkt())
    prj_file.close()

  • 相关阅读:
    oracle数据库常用指令
    MySql常用命令
    js动态添加和删除行
    Mybatis模糊查询
    laravel 成功跳转页面
    laravel 验证码
    git获取公钥和私钥以及常用的命令
    laravel css和js的引入
    git add -A 、git add -u 、 git add . 三种区别
    windows下github 出现Permission denied (publickey).解决方法
  • 原文地址:https://www.cnblogs.com/gispathfinder/p/5722103.html
Copyright © 2011-2022 走看看