命令行记录-初始Proj4包以及栅格数据投影转换

1、本篇内容包含两个部分,一是使用PROJ4包对点进行投影转换,二是栅格数据投影转换的示例

2、

#3\另外一个投影包PROJ4
from pyproj import Proj,Geod,transform

#projection1:UTM zone15, grs80 ellipse, NAD83 datum
#(defined by epsg code 26915)
p1=Proj(init=‘epsg:26915‘)
#projection2:UTM zone 15,clrk66 ellipse, NAD27 datum
p2=Proj(init=‘epsg:26715‘)
#find x,y of Jefferson City, MO.
x1,y1=p1(-92.199881,38.56694)
#transform this point to projection 2 coordinates.
x2,y2=transform(p1,p2,x1,y1)

‘%9.3f %11.3f‘ % (x1,y1)

输出:‘569704.566 4269024.671‘
‘%9.3f %11.3f‘ % (x2,y2)

输出:‘569722.342 4268814.028‘
‘%8.3f %5.3f‘ % p2(x2,y2,inverse=True)

输出:‘ -92.200 38.567‘

#process 3 points at a time in a tuple
lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri
lons = (-92.22,-94.72,-90.37)
x1, y1 = p1(lons,lats)
x2, y2 = transform(p1,p2,x1,y1)
xy=x1+y1%这里类似于配对
‘%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f‘ % xy

输出:‘567703.344 351730.944 728553.093 4298200.739 4353698.725  4292319.005‘
xy = x2+y2
‘%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f‘ % xy

输出:‘567721.149 351747.558 728569.133 4297989.112 4353489.645  4292106.305‘
lons, lats = p2(x2,y2,inverse=True)
xy = lons+lats
‘%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f‘ % xy
输出:‘ -92.220  -94.720  -90.370 38.830 39.320 38.750‘
# test datum shifting, installation of extra datum grid files.
p1 = Proj(proj=‘latlong‘,datum=‘WGS84‘)
x1 = -111.5; y1 = 45.25919444444
p2 = Proj(proj="utm",zone=10,datum=‘NAD27‘)
x2, y2 = transform(p1, p2, x1, y1)
"%s %s" % (str(x2)[:9],str(y2)[:9])
输出:‘1402291.0 5076289.5‘

#4\栅格数据投影转换
from osgeo import gdal,osr
from osgeo.gdalconst import *
#源图像投影,目标图像投影
sr1=osr.SpatialReference()
sr1.ImportFromEPSG(32650) #WGS84/UTM ZONE 50
sr2=osr.SpatialReference()
sr2.ImportFromEPSG(3857) #Google, Web Mercator
coordTrans=osr.CoordinateTransformation(sr1,sr2)

#打开源图像文件
ds1=gdal.Open("fdem.tif")
#insr=ds1.GetProjection()#WGS84/UTM ZONE 50
mat1=ds1.GetGeoTransform()
#print(mat1)
#(409294.88696681266, 27.376482012944024, 0.0, 4423871.083377095, 0.0, -27.376482012944006)

#源图像的左上角与右下角像素,在目标图像中的坐标
(ulx,uly,ulz)=coordTrans.TransformPoint(mat1[0],mat1[3])
(lrx,lry,lrz)=coordTrans.TransformPoint(mat1[0]+mat1[1]*ds1.RasterXSize,\
                                        mat1[3]+mat1[5]*ds1.RasterYSize)

#创建目标图像文件(空白图像),行列数、波段数以及数值类型仍等同原图像
driver=gdal.GetDriverByName("GTiff")
ds2=driver.Create("fdem_lonlat.tif",ds1.RasterXSize,ds1.RasterYSize,1,GDT_UInt16)

resolution=(int)((lrx-ulx)/ds1.RasterXSize)#分辨率
mat2=[ulx,resolution,0,uly,0,-resolution]#输出图像的6个仿射变换系数
ds2.SetGeoTransform(mat2)
ds2.SetProjection(sr2.ExportToWkt())#投影,文本方式

#投影转换与重采样(gdal.GRA_NearestNeighbour,gdal.GRA_Cubic,gdal.GRA_Bilinear)
gdal.ReprojectImage(ds1,ds2,sr1.ExportToWkt(),sr2.ExportToWkt(),gdal.GRA_Bilinear)

#关闭
ds1=None
ds2=None

3、投影转换结果

原文地址:https://www.cnblogs.com/vividautumn/p/11620644.html

时间: 2024-11-09 04:40:06

命令行记录-初始Proj4包以及栅格数据投影转换的相关文章

利用命令行引用外部jar包以使程序正常运行的4种方法

声明:本博客为原创博客,未经允许,不得转载!原文链接为http://blog.csdn.net/bettarwang/article/details/30976069 平时写一些小的Java Demo时我比较喜欢用UltraEdit或记事本写完后,直接利用命令行进行编译和运行.这样的好处就是方便快捷.相信有这个习惯的人应该还大有人在.但是如果要引用外部jar包,应该如何操作呢?在写JDBC的一些Demo时,由于要利用jar包来加载相应的数据库,每个Demo都用到了外部jar包,所以特地总结了一下

利用命令行引用外部jar包以使程序正常执行的4种方法

声明:本博客为原创博客.未经同意.不得转载!原文链接为http://blog.csdn.net/bettarwang/article/details/30976069 平时写一些小的Java Demo时我比較喜欢用UltraEdit或记事本写完后,直接利用命令行进行编译和执行.这种优点就是方便快捷.相信有这个习惯的人应该还大有人在. 可是假设要引用外部jar包.应该怎样操作呢?在写JDBC的一些Demo时.因为要利用jar包来载入对应的数据库.每一个Demo都用到了外部jar包,所以特地总结了一

mac命令行对复杂ipa包重新签名

最近在做ios的自动化平台,需要通过命令行安装卸载ipa包 好了问题来,别人上传的ipa包,很可能是开发签名了只能在特定手机上安装的测试ipa包,那我们如何将其安装在我们的自动化的iphone上呢? 答案看起来显而易见,将其重新签名,但是我们是自动化平台,总不能手动签名所以需要使用mac的命令行将其重新签名 这是我第一版代码: #!/bin/bash echo $1 cd uploadfiles unzip $1 rm -rf $1 cd Payload a=`ls|sed 's/[ ][ ]*

命令行参数(flag包)

命令行参数 命令行参数可以直接通过 os.Args 获取,另外标准库的 flag 包专门用于接收和解除命令行参数 os.Args 简单的只是从命令行获取一个或一组参数,可以直接使用 os.Args.下面的这种写法,无需进行判断,无论是否提供了命令行参数,或者提供了多个,都可以处理: // 把命令行参数,依次打印,每行一个 func main() { for _, s := range os.Args[1:] { fmt.Println(s) } } flag 基本使用 下面的例子使用了两种形式的

Go命令行参数解析flag包

go语言提供的flag包可以解析命令行的参数,代码: package main import ( "flag" "fmt" ) func main() { //第一个参数,为参数名称,第二个参数为默认值,第三个参数是说明 username := flag.String("name", "", "Input your username") flag.Parse() fmt.Println("Hell

命令行记录

命令行 命令行大都是英文的缩写,记住常见的英文缩写问题也就不大了 directory 目录,文件夹 file 文件 make 新建 remove 删除 move 移动 copy 复制 list 罗列 link 链接 find 查找 echo 回声,共鸣 touch 触摸 change 改变 现在将这些组合在一起 make directory 新建文件夹 --> mkdir remove 移除 rm move 移动/重命名 mv copy 复制 cp list 罗列 ls change direc

命令行编译带外部包依赖的java源文件 [以JDBC MySQL8为例]

环境: MySQL8 JDK11(SE) 首先下载MySQL8的JDBC驱动 https://dev.mysql.com/downloads/connector/j/选 PlatForm Independent 下载完文件名差不多是 “mysql-connector-java-8.<小版本号>.jar” 然后写一个JDBC小的程序(需要数据库提前建好表) 1 import java.sql.Connection; 2 import java.sql.DriverManager; 3 impor

命令行记录-python 读shp文件记录

1.读shapefile文件主要读以下内容,包括spatialRef投影信息,layerDefn图层定义信息,geomType几何对象类型,fieldDefn字段定义信息.geomlist是得到了每一个feature的geometryRef后转为Wkt形式表示的坐标点位,reclist通过name属性一次获取每个feature的所有field信息. 2. #读ArcGIS Shape文件示例from osgeo import ogrfilename="cntry98.shp" #以只读

命令行记录-csv转为shape文件

1.除了一些转为shapefile文件必要的设置外,python读取csv文件也有一些需要注意的地方. (1)首先读取了第一行fds作为字段名 (2)之后以此读取文件内容存在data里,这部分代码比较冗长 (3)代码容易出错的地方在于ds = driver.CreateDataSource(filename[:-4])这行代码,生成shapefile文件的过程,首先是创建给定名称的文件夹(我这里是stations2),然后再在文件夹里储存相应的shp等文件,文件夹的名称和子文件的名称可以不同,但