中国深圳航空公司官方网站,大朗做网站公司,枫泾网站建设,深圳创业孵化基地入驻条件1 前言
地理坐标系#xff0c;是使用三维球面来定义地球表面位置#xff0c;以实现通过经纬度对地球表面点位引用的坐标系。 地理坐标系经过地图投影操作后就变成了投影坐标系。而地图投影是按照一定的数学法则将地球椭球面上点的经维度坐标转换到平面上的直角坐标。 2 流程…1 前言
地理坐标系是使用三维球面来定义地球表面位置以实现通过经纬度对地球表面点位引用的坐标系。 地理坐标系经过地图投影操作后就变成了投影坐标系。而地图投影是按照一定的数学法则将地球椭球面上点的经维度坐标转换到平面上的直角坐标。 2 流程
2.1 矢量数据地理坐标转投影坐标
要素feature的形状geometry是由一系列点坐标构成的将每个要素的形状点一一进行坐标转换即可。
下面以WGS84坐标转UTM投影为例
from osgeo import ogr,osr,gdal
import glob
import osdef vecter_WGS2UTM(shp_path, UTM_shp_path, longitude, is_north):ds ogr.Open(shp_path)layer ds.GetLayer(0)driver ogr.GetDriverByName(ESRI Shapefile)# 创建输出文件if os.path.exists(UTM_shp_path):driver.DeleteDataSource(UTM_shp_path)out_ds driver.CreateDataSource(UTM_shp_path)outlayer out_ds.CreateLayer(UTM_shp_path[:-4],geom_type ogr.wkbPolygon)# 当前地理参考spatialRef layer.GetFeature(0).GetGeometryRef().GetSpatialReference()# 根据经度计算UTM区号,进而定义UTM投影zone str(int(longitude/6) 31)zone int(326 zone) if is_north else int(327 zone)UTM_spatialRef osr.SpatialReference()UTM_spatialRef.ImportFromEPSG(zone)# 投影转换coordinate_transfor osr.CoordinateTransformation(spatialRef, UTM_spatialRef)# 定义输出属性表信息feature layer.GetFeature(0)field_count feature.GetFieldCount()field_names []for attr in range(field_count):field_defn feature.GetFieldDefnRef(attr)field_names.append(field_defn.GetName())outlayer.CreateField(field_defn)out_fielddefn outlayer.GetLayerDefn()for feature in layer:geometry feature.GetGeometryRef()geometry.Transform(coordinate_transfor)out_feature ogr.Feature(out_fielddefn)out_feature.SetGeometry(geometry)for field_name in field_names:out_feature.SetField(field_name,feature.GetField(field_name))outlayer.CreateFeature(out_feature)feature.Destroy()out_feature.Destroy()# 清除缓存ds.Destroy()out_ds.Destroy()# 保存投影文件UTM_spatialRef.MorphFromESRI()prj_path UTM_shp_path.replace(.shp,.prj)fn open(prj_path,w)fn.write(UTM_spatialRef.ExportToWkt())fn.close()
2.2 栅格数据地理坐标转投影坐标
栅格数据每个像素的地理/投影坐标是由仿射矩阵六参数和像素坐标计算得来的所以先将仿射矩阵六参数进行转换之后对栅格数据重采样即可。
下面以WGS84坐标转UTM投影为例
def raster_WGS2UTM(raster_path, UTM_raster_path, longitude, is_north):raster_ds gdal.Open(raster_path)raster_type raster_ds.GetRasterBand(1).DataType# 栅格投影spatialRef osr.SpatialReference()spatialRef.ImportFromWkt(raster_ds.GetProjection())# 根据经度计算UTM区号,进而定义UTM投影zone str(int(longitude/6) 31)zone int(326 zone) if is_north else int(327 zone)UTM_spatialRef osr.SpatialReference()UTM_spatialRef.ImportFromEPSG(zone)# 投影转换coordinate_transfor osr.CoordinateTransformation(spatialRef, UTM_spatialRef)# 仿射矩阵六参数geotransform raster_ds.GetGeoTransform()# 左上角upper left、右下角lower right坐标ul_x geotransform[0]ul_y geotransform[3]lr_x geotransform[0]geotransform[1]*raster_ds.RasterXSizegeotransform[2]*raster_ds.RasterYSizelr_y geotransform[3]geotransform[4]*raster_ds.RasterYSizegeotransform[5]*raster_ds.RasterYSize# 左上角、右下角在目标投影中的坐标(UTM_ul_x, UTM_ul_y, UTM_ul_z) coordinate_transfor.TransformPoint(ul_y, ul_x)(UTM_lr_x, UTM_lr_y, UTM_lr_z) coordinate_transfor.TransformPoint(lr_y, lr_x)# 创建目标图像文件driver gdal.GetDriverByName(GTiff)UTM_raster_ds driver.Create(UTM_raster_path, raster_ds.RasterXSize,raster_ds.RasterYSize,raster_ds.RasterCount,raster_type)# 转换后图像的分辨率resolution (UTM_lr_x-UTM_ul_x)/raster_ds.RasterXSize# 转换后图像的六个放射变换参数UTM_transform [UTM_ul_x, resolution, 0, UTM_ul_y, 0, -resolution]UTM_raster_ds.SetGeoTransform(UTM_transform)UTM_raster_ds.SetProjection(UTM_spatialRef.ExportToWkt())# 投影转换后需要做重采样gdal.ReprojectImage(raster_ds, UTM_raster_ds, spatialRef.ExportToWkt(), UTM_spatialRef.ExportToWkt(), gdal.GRA_Bilinear)# 关闭raster_ds NoneUTM_raster_ds None
来源应用推广部
供稿技术研发部
编辑方梅
