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()

  • 相关阅读:
    高效App运营必须get的七个推送技巧
    如何预测您的用户即将流失?
    HMS Core赋能移动金融,为行业注入增长新动力
    做好自定义预测,探寻产品增长动能
    Unity平台 | 快速集成华为AGC云数据库服务
    HMS Core Insights第一期直播回顾 – 深入浅出,创新技术与开发者共同推进行业发展!
    多种细分方式浏览销售数据,IAP助您有效洞察市场收益效果
    Unity | 快速集成华为AGC云存储服务
    awk命令
    反转字符串
  • 原文地址:https://www.cnblogs.com/gispathfinder/p/5722103.html
Copyright © 2011-2022 走看看