有时候需要获取遥感影像的范围,即落图文件。利用高分影像的元数据文件生成则十分简便快捷。根据实际需求,编写了生成落图的程序,主要针对高分的影像。不同的传感器的数据存储的形式也不一定一样,针对相同的存储形式写了相应的函数。基于Python实现。考虑到不同的影像有不同的投影,最后每个xml文件对应一个shp。
1 导入需要的库
import os
import osgeo.ogr as ogr
import osgeo.osr as osr
from xml.etree import ElementTree
import tkFileDialog
2 设置输入输出路径 ,手动设置或用导航式选择
#inputpath="E:\\wcg\\luotu"
#outputpath="E:\\wcg\\luotu13"
inputpath=tkFileDialog.askdirectory(title="Select inputs XML folder")
outputpath=tkFileDialog.askdirectory(title="Select outputs SHPFILE folder")
3 获取输入路径内的所有文件
def getListFiles(path):
assert os.path.isdir(path),'%s not exist,'%path
ret=[]
for root,dirs,files in os.walk(path):
for filespath in files:
ret.append(os.path.join(root,filespath))
return ret
4 针对ZY3、02C、GF2、GF1、SV1 的函数
def read_ZY3_02C_GF2_GF1_SV1(text):
corner={}
corner["A"]=[]
corner["B"]=[]
corner["C"]=[]
corner["D"]=[]
others={}
others["SolarAzimuth"]=""
others["SolarZenith"]=""
others["SatelliteAzimuth"]=""
others["SatelliteZenith"]=""
others["StartTime"]=""
others["EndTime"]=""
others["CloudPercent"]=""
root = ElementTree.parse(text)
lst_node = root.getiterator("TopLeftLatitude")
for node in lst_node:
corner["A"].append(node.text)
lst_node = root.getiterator("TopLeftLongitude")
for node in lst_node:
corner["A"].append(node.text)
lst_node = root.getiterator("TopRightLatitude")
for node in lst_node:
corner["B"].append(node.text)
lst_node = root.getiterator("TopRightLongitude")
for node in lst_node:
corner["B"].append(node.text)
lst_node = root.getiterator("BottomRightLatitude")
for node in lst_node:
corner["C"].append(node.text)
lst_node = root.getiterator("BottomRightLongitude")
for node in lst_node:
corner["C"].append(node.text)
lst_node = root.getiterator("BottomLeftLatitude")
for node in lst_node:
corner["D"].append(node.text)
lst_node = root.getiterator("BottomLeftLongitude")
for node in lst_node:
corner["D"].append(node.text)
lst_node = root.getiterator("SolarAzimuth")
for node in lst_node:
others["SolarAzimuth"]=node.text
lst_node = root.getiterator("SolarZenith")
for node in lst_node:
others["SolarZenith"]=node.text
lst_node = root.getiterator("SatelliteAzimuth")
for node in lst_node:
others["SatelliteAzimuth"]=node.text
lst_node = root.getiterator("SatelliteZenith")
for node in lst_node:
others["SatelliteZenith"]=node.text
lst_node = root.getiterator("StartTime")
for node in lst_node:
others["StartTime"]=node.text
lst_node = root.getiterator("EndTime")
for node in lst_node:
others["EndTime"]=node.text
lst_node = root.getiterator("CloudPercent")
for node in lst_node:
others["CloudPercent"]=node.text
return corner,others
5 针对YG24的函数
def read_YG24(text):
per=ElementTree.parse(text)
corners=per.findall('Corners')
four={}
four["A"]=[]
four["B"]=[]
four["C"]=[]
four["D"]=[]
others={}
others["SolarAzimuth"]=""
others["SolarZenith"]=""
others["SatelliteAzimuth"]=""
others["SatelliteZenith"]=""
others["StartTime"]=""
others["EndTime"]=""
others["CloudPercent"]=""
for corner in corners:
for child in corner.getchildren():
if child.tag=="UpperLeft":
for chi in child:
four["A"].append(chi.text)
if child.tag=="UpperRight":
for chi in child:
four["B"].append(chi.text)
if child.tag=="LowerRight":
for chi in child:
four["C"].append(chi.text)
if child.tag=="LowerLeft":
for chi in child:
four["D"].append(chi.text)
metas=per.findall('ProductMeta')
for meta in metas:
for child in meta.getchildren():
if child.tag=="StartTime":
for chi in child:
others["StartTime"]=chi.text
if child.tag=="EndTime":
for chi in child:
others["EndTime"]=chi.text
return four,others
6 针对SPOT6的函数
def read_SPOT6(text):
per=ElementTree.parse(text)
content=per.findall('Dataset_Content')
extent=content[0].findall("Dataset_Extent")
fourlist=[]
for corner in extent:
for child in corner.getchildren():
if child.tag=="Vertex":
for chi in child:
if chi.tag=="LON":
fourlist.append(chi.text)
if chi.tag=="LAT":
fourlist.append(chi.text)
four={}
four["A"]=[fourlist[1],fourlist[0]]
four["B"]=[fourlist[3],fourlist[2]]
four["C"]=[fourlist[5],fourlist[4]]
four["D"]=[fourlist[7],fourlist[6]]
others={}
others["SolarAzimuth"]=""
others["SolarZenith"]=""
others["SatelliteAzimuth"]=""
others["SatelliteZenith"]=""
others["StartTime"]=""
others["EndTime"]=""
others["CloudPercent"]=""
cloud=(content[0].findall("CLOUD_COVERAGE"))[0].text
others["CloudPercent"]=cloud
sources=per.findall('Dataset_Sources')
times=sources[0].findall("Source_Identification")
for time in times:
for child in time.getchildren():
if child.tag=="Strip_Source":
for chi in child:
if chi.tag=="IMAGING_DATE":
others["StartTime"]=chi.text+" "
if chi.tag=="IMAGING_TIME":
others["StartTime"]+=chi.text
others["EndTime"]= others["StartTime"]
return four,others
7 转出到矢量的函数
def toshp(ABCD1,otherss,outShapefile):
ring = ogr.Geometry(ogr.wkbLinearRing)
for jw in ABCD1:
ring.AddPoint(jw[0],jw[1])
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
outDriver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(outShapefile):
outDriver.DeleteDataSource(outShapefile)
outDataSource = outDriver.CreateDataSource(outShapefile)
outLayer = outDataSource.CreateLayer("states_extent",srs,geom_type=ogr.wkbPolygon)
yxmc=ogr.FieldDefn("YXMC", ogr.OFTString)
yxmc.SetWidth(255)
outLayer.CreateField(ogr.FieldDefn("ID", ogr.OFTInteger))
outLayer.CreateField(yxmc)
SolarAzimuth=ogr.FieldDefn("SolAzimuth", ogr.OFTString)
SolarAzimuth.SetWidth(10)
outLayer.CreateField(SolarAzimuth)
SolarZenith=ogr.FieldDefn("SolZenith", ogr.OFTString)
SolarZenith.SetWidth(10)
outLayer.CreateField(SolarZenith)
SatelliteAzimuth=ogr.FieldDefn("SatAzimuth", ogr.OFTString)
SatelliteAzimuth.SetWidth(10)
outLayer.CreateField(SatelliteAzimuth)
SatelliteZenith=ogr.FieldDefn("SatZenith", ogr.OFTString)
SatelliteZenith.SetWidth(10)
outLayer.CreateField(SatelliteZenith)
StartTime=ogr.FieldDefn("StartTime", ogr.OFTString)
StartTime.SetWidth(30)
outLayer.CreateField(StartTime)
EndTime=ogr.FieldDefn("EndTime", ogr.OFTString)
EndTime.SetWidth(30)
outLayer.CreateField(EndTime)
CloudPercent=ogr.FieldDefn("CloudPer", ogr.OFTString)
CloudPercent.SetWidth(10)
outLayer.CreateField(CloudPercent)
featureDefn = outLayer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(poly)
feature.SetField("id", 1)
index1=outShapefile.rfind("/")
feature.SetField("YXMC", outShapefile[index1+1:-4])
feature.SetField("SolAzimuth",otherss["SolarAzimuth"])
feature.SetField("SolZenith",otherss["SolarZenith"])
feature.SetField("SatAzimuth",otherss["SatelliteAzimuth"])
feature.SetField("SatZenith",otherss["SatelliteZenith"])
feature.SetField("StartTime",otherss["StartTime"])
feature.SetField("EndTime",otherss["EndTime"])
feature.SetField("CloudPer",otherss["CloudPercent"])
outLayer.CreateFeature(feature)
feature = None
outDataSource = None
8 最后的调用
ListFiles=getListFiles(inputpath)
b=[]
for i in ListFiles:
if (i[-3:].lower()=="xml") and ("order" not in i) and ("shp" not in i)\
and ("aux" not in i)and ("tiff" not in i)and ("orientation" not in i) :
b.append(i)
for f in b:
if f.find("DIM_SPOT6") >= 0 :
ABCD,otherss=read_SPOT6(f)
index=f.find("DIM_SPOT6")
elif f.find("DIM_SPOT7") >= 0 :
ABCD,otherss=read_SPOT6(f)
index=f.find("DIM_SPOT7")
elif f.find("YG24") >= 0 :
ABCD,otherss=read_YG24(f)
index=f.find("YG24")
elif (f.find("GF1") >=0):
ABCD,otherss=read_ZY3_02C_GF2_GF1_SV1(f)
index=f.find("GF1")
elif f.find("GF2") >=0:
ABCD,otherss=read_ZY3_02C_GF2_GF1_SV1(f)
index=f.find("GF2")
elif f.find("ZY02C") >=0:
ABCD,otherss=read_ZY3_02C_GF2_GF1_SV1(f)
index=f.find("ZY02C")
elif f.find("ZY3") >=0 :
ABCD,otherss=read_ZY3_02C_GF2_GF1_SV1(f)
index=f.find("ZY3")
elif f.find("SV1") >=0 :
ABCD,otherss=read_ZY3_02C_GF2_GF1_SV1(f)
index=f.find("SV1")
else:
continue
for i in ABCD.keys():
ABCD[i].reverse()
ABCD1=[ABCD["A"],ABCD["B"],ABCD["C"],ABCD["D"],ABCD["A"]]
for jingwei in ABCD1:
for j in range(len(jingwei)):
jingwei[j]=float(jingwei[j])
name=f[index:-4]
if not os.path.exists(outputpath+"\\"+name+".shp"):
toshp(ABCD1,otherss,outputpath+"\\"+name+".shp")
利用xml生成落图的到此结束。
下面是生成tiff文件落图的,大致思路同上,只不过获取四个角点的坐标由xml读取变为从tfw文件读取。读tfw与读普通的文本文件一致的方法,只针对tif格式影像。
扫描二维码关注公众号,回复:
2257443 查看本文章
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 26 16:23:14 2018
北京落图
https://desktop.arcgis.com/zh-cn/arcmap/latest/manage-data/raster-and-images/world-files-for-raster-datasets.htm
@author: Administrator
"""
import os
import osgeo.ogr as ogr
import osgeo.osr as osr
import gdal
import tkFileDialog
#inputpath=tkFileDialog.askdirectory(title="input")
#outputpath=tkFileDialog.askdirectory(title="output")
inputpath="E:\\hb哈哈哈哈".decode("utf-8")
outputpath="E:\\wcg\\aa"
target=outputpath
def getListFiles(path):
assert os.path.isdir(path),'%s not exist,'%path
ret=[]
for root,dirs,files in os.walk(path):
for filespath in files:
ret.append(os.path.join(root,filespath))
return ret
inputpat=inputpath.replace("/","\\")
files=getListFiles(inputpat)
def image(path):
fi=path[:-3]+"TFW"
f=open(fi,"r")
c=f.read().split("\n")
f.close()
A,D,B,E,C,F=float(c[0]),float(c[1]),float(c[2]),float(c[3]),float(c[4]),float(c[5])
dataset = gdal.Open(path)
nXSize = dataset.RasterXSize
nYSize = dataset.RasterYSize
geotransform = dataset.GetGeoTransform()
dx=geotransform[1]
dy=geotransform[5]
projinfo=dataset.GetProjection()
Ax,Ay=A*-0.5+B*-0.5+C,D*-0.5+E*-0.5+F
Bx,By=A*(nXSize-0.5)+B*(-0.5)+C,D*(nXSize-0.5)+E*-0.5+F
Cx,Cy=A*(nXSize-0.5)+B*(nYSize-0.5)+C,D*(nXSize-0.5)+E*(nYSize-0.5)+F
Dx,Dy=A*-0.5+B*(nYSize-0.5)+C,D*-0.5+E*(nYSize-0.5)+F
four={}
four["A"]=[Ax,Ay]
four["B"]=[Bx,By]
four["C"]=[Cx,Cy]
four["D"]=[Dx,Dy]
return four,projinfo
def toshp(ABCD1,outShapefile,projinfo,lujing):
ring = ogr.Geometry(ogr.wkbLinearRing)
for jw in ABCD1:
ring.AddPoint(jw[0],jw[1])
srs = osr.SpatialReference()
srs.SetFromUserInput( projinfo )
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
outDriver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(outShapefile):
outDriver.DeleteDataSource(outShapefile)
outDataSource = outDriver.CreateDataSource(outShapefile)
outLayer = outDataSource.CreateLayer("states_extent",srs,geom_type=ogr.wkbPolygon)
yxmc=ogr.FieldDefn("YXMC", ogr.OFTString)
yxmc.SetWidth(255)
outLayer.CreateField(ogr.FieldDefn("ID", ogr.OFTInteger))
outLayer.CreateField(yxmc)
lj=ogr.FieldDefn("lujing", ogr.OFTString)
lj.SetWidth(255)
outLayer.CreateField(lj)
featureDefn = outLayer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(poly)
feature.SetField("id", 1)
index1=outShapefile.rfind("\\")
feature.SetField("YXMC", outShapefile[index1+1:-4])
feature.SetField("lujing", lujing)
outLayer.CreateFeature(feature)
feature = None
outDataSource = None
jpfs=[]
#os.mkdir(target)
for i in files:
if i[-3:].lower() in ["tif","tiff"]:
jpfs.append(i)
for a in jpfs:
ABCD=image(a)[0]
ABCD1=[ABCD["A"],ABCD["B"],ABCD["C"],ABCD["D"],ABCD["A"]]
for jingwei in ABCD1:
for j in range(len(jingwei)):
jingwei[j]=float(jingwei[j])
index1=a.rfind("\\")
name=a[index1:-4]
lujing=a[:index1+1]
toshp(ABCD1,target+"\\"+name+".shp",image(a)[1],lujing)