使用 Python 批量获取遥感影像的质心点坐标并对应影像

(x + RasterXSize * x_pixel_resolution/2 , y +RasterYSize * y_pixel_resolution/2)

复制

from osgeo import gdal
import os
import csv

#initialize
image_dir = './my_dir/'
image_files = os.listdir(image_dir)
image_files.sort(key=lambda x:int(x[:-4]))

with open('test.csv','w',encoding='utf-8',newline='') as f:
    writer = csv.writer(f)
    writer.writerow(["file_name","x_mid","y_mid"])
#读取影像
    for file in image_files:

        image_file = os.path.join(image_dir,file)

        in_ds = gdal.Open(image_file)
        in_ds_trans = in_ds.GetGeoTransform()

#宽度和分辨率
        width = in_ds.RasterXSize
        height = in_ds.RasterYSize
        width_pixel = in_ds_trans[1]
        height_pixel = in_ds_trans[5]

#计算中心点的坐标
        x_left_top = in_ds_trans[0]
        y_left_top = in_ds_trans[3]

        x_mid = x_left_top + width/2 * width_pixel
        y_mid = y_left_top + height/2 * height_pixel

        row_info = [file,x_mid,y_mid]
        writer.writerow(row_info)
        del in_ds

复制

import pandas as pd
import numpy as np
import csv
f = open('find_file.csv','w',encoding='utf-8',newline='')
csv_writer = csv.writer(f)
data = pd.read_csv('./随便取名就好了.csv',header=None)
data1 = pd.read_csv('./取名开心就好了.csv',header = None)

#省略表头 读取 attribute
net_x = data.iloc[:,2][1:len(data)]
net_y = data.iloc[:,3][1:len(data)]
completeness = data.iloc[:,6][1:len(data)]

Name = data1.iloc[:,0][1:len(data1)]
x = data1.iloc[:,1][1:len(data1)]
y = data1.iloc[:,2][1:len(data1)]

#x_mid_data 表示的是要找的点
#转 list 并且转成 float 类型
net_x = list(net_x)
net_x_data = [float(data) for data in net_x]
net_y = list(net_y)
net_y_data = [float(data) for data in net_y]
completeness = list(completeness)
completeness_data = [float(data) for data in atttibute]

# 要找的图片
Name = list(Name)
x = list(x)
y = list(y)

x_data = [float(data) for data in x]
y_data = [float(data) for data in y]


#  代表的是 test.csv,图片计算出来的点
# x_data[i] 代表的是计算出来的质心点

for i in range(len(net_x)):
    for index,xy in enumerate(zip(x_data,y_data)):
        if  abs(xy[0] -net_x_data[i])<1  and abs(xy[1] - net_y_data[i])<1:
            write_info = [xy[0],xy[1],Name[index],attribute[i]]
            csv_writer.writerow(write_info)

&emsp;&emsp;思路其实很简单，就是用 ArcGIS 的生成的表里面的质心去搜索另外一张表，如果找到的话就返回它的 `index`,所以外层循环就是 ArcGIS 的表里面搜索出的质心(x,y),然后在另一个表内搜索相应的位置，返回坐标。


```python
a = [1,2,3]
b = [2,3,4]
c = zip(a,b)
print(list(c))
d = [1,2,3]
e = [2,3,4]

for i in range(3):
    for index,xy in enumerate(zip(a,b)):
        if a[i] == xy[0] and b[i] == xy[1]:
            print(f'坐标对应的索引是{index}')
```

output:
    [(1, 2), (2, 3), (3, 4)]
    坐标对应的索引是0
    坐标对应的索引是1
    坐标对应的索引是2

复制