Skip to content

Commit 785646e

Browse files
ZivKiddZivKidd
committedNov 22, 2020
first commit
1 parent 35b50ae commit 785646e

28 files changed

+4650
-0
lines changed
 

‎1duxiang.py

+93
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
# -*- coding: utf-8 -*-
2+
import gdal
3+
import numpy as np
4+
cell=1024
5+
clip_folder=r"H:\xzr\duxiang_buffer\clip_1024\\"
6+
7+
# 读取要切的原图
8+
in_ds = gdal.Open(r"H:\xzr\duxiang_buffer\2kmclip.tif")
9+
print("open tif file succeed")
10+
11+
# 读取原图中的每个波段
12+
in_band1 = in_ds.GetRasterBand(1)
13+
in_band2 = in_ds.GetRasterBand(2)
14+
in_band3 = in_ds.GetRasterBand(3)
15+
16+
w=in_ds.RasterXSize
17+
h=in_ds.RasterYSize
18+
19+
w=int(w/(cell/2))-1
20+
h=int(h/(cell/2))-1
21+
22+
for i in range(w):
23+
for j in range(h):
24+
# 定义切图的起始点坐标(相比原点的横坐标和纵坐标偏移量)
25+
offset_x =int(i*cell/2)
26+
offset_y =int(j*cell/2)
27+
28+
29+
# 从每个波段中切需要的矩形框内的数据(注意读取的矩形框不能超过原图大小)
30+
out_band1 = in_band1.ReadAsArray(offset_x, offset_y, cell, cell)
31+
out_band2 = in_band2.ReadAsArray(offset_x, offset_y, cell, cell)
32+
out_band3 = in_band3.ReadAsArray(offset_x, offset_y, cell, cell)
33+
34+
if(np.sum(out_band1)==0 and np.sum(out_band2)==0 and np.sum(out_band3)==0):
35+
continue
36+
37+
if(np.where(out_band1==0)[0].shape[0]>1024**2/2):
38+
continue
39+
40+
print(offset_x,offset_y)
41+
42+
43+
# 获取Tif的驱动，为创建切出来的图文件做准备
44+
gtif_driver = gdal.GetDriverByName("GTiff")
45+
46+
# 创建切出来的要存的文件（3代表3个不都按，最后一个参数为数据类型，跟原文件一致）
47+
out_ds = gtif_driver.Create(clip_folder+str(offset_x)+'-'+str(offset_y)+ '.tif', cell, cell, 3, in_band1.DataType)
48+
print("create new tif file succeed")
49+
50+
# 获取原图的原点坐标信息
51+
ori_transform = in_ds.GetGeoTransform()
52+
if ori_transform:
53+
print (ori_transform)
54+
print("Origin = ({}, {})".format(ori_transform[0], ori_transform[3]))
55+
print("Pixel Size = ({}, {})".format(ori_transform[1], ori_transform[5]))
56+
57+
# 读取原图仿射变换参数值
58+
top_left_x = ori_transform[0] # 左上角x坐标
59+
w_e_pixel_resolution = ori_transform[1] # 东西方向像素分辨率
60+
top_left_y = ori_transform[3] # 左上角y坐标
61+
n_s_pixel_resolution = ori_transform[5] # 南北方向像素分辨率
62+
63+
# 根据反射变换参数计算新图的原点坐标
64+
top_left_x = top_left_x + offset_x * w_e_pixel_resolution
65+
top_left_y = top_left_y + offset_y * n_s_pixel_resolution
66+
67+
# 将计算后的值组装为一个元组，以方便设置
68+
dst_transform = (top_left_x, ori_transform[1], ori_transform[2], top_left_y, ori_transform[4], ori_transform[5])
69+
70+
# 设置裁剪出来图的原点坐标
71+
out_ds.SetGeoTransform(dst_transform)
72+
73+
# 设置SRS属性（投影信息）
74+
out_ds.SetProjection(in_ds.GetProjection())
75+
76+
# 写入目标文件
77+
out_ds.GetRasterBand(1).WriteArray(out_band1)
78+
out_ds.GetRasterBand(2).WriteArray(out_band2)
79+
out_ds.GetRasterBand(3).WriteArray(out_band3)
80+
81+
# 将缓存写入磁盘
82+
out_ds.FlushCache()
83+
print("FlushCache succeed")
84+
85+
# 计算统计值
86+
# for i in range(1, 3):
87+
# out_ds.GetRasterBand(i).ComputeStatistics(False)
88+
# print("ComputeStatistics succeed")
89+
90+
del out_ds
91+
92+
print("End!")
93+

‎1read_tif.py

+169
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,169 @@
1+
# -*- coding: utf-8 -*-
2+
import gdal
3+
import osgeo
4+
import os
5+
import shutil
6+
import shapefile
7+
from osgeo import osr
8+
import numpy as np
9+
import arcpy
10+
from PIL import Image, ImageDraw
11+
12+
arcpy.env.workspace = r"D:\xzr\process\data.gbd" # arcgis地理数据库目录
13+
14+
15+
def lonlat2geo(dataset, lon, lat):
16+
'''
17+
将经纬度坐标转为投影坐标（具体的投影坐标系由给定数据确定）
18+
:param dataset: GDAL地理数据
19+
:param lon: 地理坐标lon经度
20+
:param lat: 地理坐标lat纬度
21+
:return: 经纬度坐标(lon, lat)对应的投影坐标
22+
'''
23+
prosrs, geosrs = getSRSPair(dataset)
24+
ct = osr.CoordinateTransformation(geosrs, prosrs)
25+
coords = ct.TransformPoint(lon, lat)
26+
return coords[:2]
27+
28+
29+
def getSRSPair(dataset):
30+
'''
31+
获得给定数据的投影参考系和地理参考系
32+
:param dataset: GDAL地理数据
33+
:return: 投影参考系和地理参考系
34+
'''
35+
prosrs = osr.SpatialReference()
36+
prosrs.ImportFromWkt(dataset.GetProjection())
37+
geosrs = prosrs.CloneGeogCS()
38+
return prosrs, geosrs
39+
40+
41+
def geo2imagexy(dataset, x, y):
42+
'''
43+
根据GDAL的六 参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
44+
:param dataset: GDAL地理数据
45+
:param x: 投影或地理坐标x
46+
:param y: 投影或地理坐标y
47+
:return: 影坐标或地理坐标(x, y)对应的影像图上行列号(row, col)
48+
'''
49+
trans = dataset.GetGeoTransform()
50+
a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
51+
b = np.array([x - trans[0], y - trans[3]])
52+
return np.linalg.solve(a, b) # 使用numpy的linalg.solve进行二元一次方程的求解
53+
54+
55+
def lonlat2imagexy(dataset, x, y):
56+
'''
57+
影像行列转经纬度：
58+
：通过经纬度转平面坐标
59+
：平面坐标转影像行列
60+
'''
61+
coords = lonlat2geo(dataset, x, y)
62+
coords2 = geo2imagexy(dataset, coords[0], coords[1])
63+
return (int(round(abs(coords2[0]))), int(round(abs(coords2[1]))))
64+
65+
66+
def getAllFileName(folder_path):
67+
file_list = []
68+
folder_list = []
69+
for file_name in os.listdir(folder_path):
70+
if (os.path.isfile(os.path.join(folder_path, file_name))):
71+
file_list.append(os.path.join(folder_path, file_name))
72+
elif (os.path.isdir(os.path.join(folder_path, file_name))):
73+
folder_list.append(os.path.join(folder_path, file_name))
74+
file_list.sort()
75+
return file_list, folder_list
76+
77+
78+
def getFileName(file_dir):
79+
file_path_list = []
80+
for root, dirs, files in os.walk(file_dir):
81+
for file in files:
82+
if (file[-3:] == "tif"):
83+
file_path_list.append(unicode(root + '\\' + file,'gbk'))
84+
return file_path_list
85+
86+
87+
tif_folders = [r"D:\deepleearning\sampleshp\tif_all"]
88+
tif_files = []
89+
for t in tif_folders:
90+
tif_files.extend(getFileName(t))
91+
92+
shp_path = r"D:\deepleearning\sampleshp\titian\titian.shp" # shp文件的路径， shapefile不支持中文路径
93+
out_dir = r"D:\deepleearning\sampleshp\titian\tif_clip" # 裁剪后图像保存路径
94+
sf = shapefile.Reader(shp_path) # 读取shp文件
95+
shapes = sf.shapes()
96+
97+
for i in range(len(shapes)):
98+
print str(i) + '/' + str(len(shapes))
99+
shp = shapes[i] # 获取shp文件中的每一个形状
100+
101+
point = shp.points # 获取每一个最小外接矩形的四个点
102+
x_list = [ii[0] for ii in point]
103+
y_list = [ii[1] for ii in point]
104+
105+
x_min = min(x_list)
106+
y_min = min(y_list)
107+
x_max = max(x_list)
108+
y_max = max(y_list)
109+
110+
x_cen = (x_min + x_max) / 2
111+
y_cen = (y_max + y_min) / 2
112+
113+
x_min1 = x_min - (x_max - x_min) / 2
114+
x_max1 = x_max + (x_max - x_min) / 2
115+
y_min1 = y_min - (y_max - y_min) / 2
116+
y_max1 = y_max + (y_max - y_min) / 2
117+
#
118+
# x_min1 = x_min - (x_max - x_min) * 2
119+
# x_max1 = x_max + (x_max - x_min) * 2
120+
# y_min1 = y_min - (y_max - y_min) * 2
121+
# y_max1 = y_max + (y_max - y_min) * 2
122+
123+
# x_min1 = x_min - (x_max - x_min)
124+
# x_max1 = x_max + (x_max - x_min)
125+
# y_min1 = y_min - (y_max - y_min)
126+
# y_max1 = y_max + (y_max - y_min)
127+
128+
count = -1
129+
for t in tif_files:
130+
dataset = gdal.Open(t)
131+
im_width = dataset.RasterXSize # 栅格矩阵的列数
132+
im_height = dataset.RasterYSize # 栅格矩阵的行数
133+
134+
im_geotrans = dataset.GetGeoTransform() # 仿射矩阵
135+
im_proj = dataset.GetProjection() # 地图投影信息
136+
137+
im_x_min = im_geotrans[0]
138+
im_y_max = im_geotrans[3]
139+
im_x_max = im_x_min + im_width * im_geotrans[1]
140+
im_y_min = im_y_max + im_height * im_geotrans[5]
141+
142+
if (y_cen < im_y_max and y_cen > im_y_min and x_cen < im_x_max and x_cen > im_x_min):
143+
144+
coords = lonlat2imagexy(dataset, x_cen, y_cen)
145+
coords1 = lonlat2geo(dataset, x_max1, y_max1)
146+
coords2 = lonlat2geo(dataset, x_min1, y_min1)
147+
x_min2 = coords2[0]
148+
y_min2 = coords2[1]
149+
x_max2 = coords1[0]
150+
y_max2 = coords1[1]
151+
152+
if (coords[0] < im_width and coords[1] < im_height):
153+
print t
154+
count += 1
155+
156+
out_path = os.path.join(out_dir, str(i) + '-' + str(count) + '.tif')
157+
158+
cor = str(x_min2) + ' ' + str(y_min2) + ' ' + str(x_max2) + ' ' + str(y_max2)
159+
# cor = str(x_min1) + ' ' + str(y_min1) + ' ' + str(x_max1) + ' ' + str(y_max1)
160+
print cor
161+
arcpy.Clip_management(t, cor, out_path, "#", "#", None) # 调用工具箱函数
162+
163+
try:
164+
img = Image.open(out_path)
165+
print 'good'
166+
except:
167+
os.remove(out_path)
168+
os.remove(out_path + '.ovr')
169+
print 'deleted'

‎2duxiang.py

+71
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
import gdal
2+
import numpy as np
3+
import os
4+
from PIL import Image, ImageDraw
5+
import datetime
6+
7+
def getAllFileName(folder_path):
8+
file_list = []
9+
folder_list = []
10+
for file_name in os.listdir(folder_path):
11+
if (os.path.isfile(os.path.join(folder_path, file_name))):
12+
file_list.append(os.path.join(folder_path, file_name))
13+
elif (os.path.isdir(os.path.join(folder_path, file_name))):
14+
folder_list.append(os.path.join(folder_path, file_name))
15+
file_list.sort()
16+
return file_list, folder_list
17+
18+
19+
# raster = np.zeros((1000,2000), dtype=np.int)
20+
# raster[200:300,1000:1800]=255
21+
#
22+
# # 读取要切的原图
23+
# in_ds = gdal.Open(r"H:\xzr\duxiang\压缩文件\201912_GF2.img")
24+
# in_band1 = in_ds.GetRasterBand(1)
25+
#
26+
# dst_ds = gdal.GetDriverByName('GTiff').Create("hello.tif", 2000, 1000, 1, in_band1.DataType)
27+
# dst_ds.SetGeoTransform(in_ds.GetGeoTransform())
28+
# dst_ds.SetProjection(in_ds.GetProjection())
29+
# dst_ds.GetRasterBand(1).WriteArray(raster)
30+
# # Once we're done, close properly the dataset
31+
# dst_ds.FlushCache()
32+
33+
34+
35+
raster = np.zeros((78234,147474), dtype=np.int)
36+
37+
38+
file_list,_=getAllFileName(r"D:\semantic-segmentation-for-Geographical-survey\Val1")
39+
num=0
40+
for f in file_list:
41+
num += 1
42+
if(num%1000==0):
43+
print(datetime.datetime.now())
44+
print(num,'/25731')
45+
46+
name=os.path.split(f)[1][:-4]
47+
name=name.split('-')
48+
x=int(name[0])
49+
y=int(name[1])
50+
51+
image=Image.open(f)
52+
image=image.resize((1024,1024))
53+
image = np.array(image) # image类 转 numpy
54+
image = image[:, :, 0] # 第1通道
55+
56+
if(np.max(image)==0):
57+
continue
58+
59+
has_value=np.where(image>0)
60+
for v in range(has_value[0].shape[0]):
61+
raster[has_value[0][v]+y,has_value[1][v]+x]=255
62+
63+
in_ds = gdal.Open(r"H:\xzr\duxiang\压缩文件\201912_GF2.img")
64+
in_band1 = in_ds.GetRasterBand(1)
65+
66+
dst_ds = gdal.GetDriverByName('GTiff').Create("hello4.tif", 147474, 78234, 1, in_band1.DataType)
67+
dst_ds.SetGeoTransform(in_ds.GetGeoTransform())
68+
dst_ds.SetProjection(in_ds.GetProjection())
69+
dst_ds.GetRasterBand(1).WriteArray(raster)
70+
# Once we're done, close properly the dataset
71+
dst_ds.FlushCache()

0 commit comments

Comments
 (0)
Please sign in to comment.