python实现对遥感影像经纬度获取并实现海陆分离

1 珊格数据（.pix、.lyr等）转tif

tif百科
转tif

2 使用gdal库实现经纬度获取

tif图像自带经纬度信息，因此可以使用gdal包获取经纬度信息然后对相应像素点负值操作

from osgeo import gdal
from global_land_mask import globe
from  global_land_mask.globe import is_land
from glob import *
from osgeo import osr
import os
# import numpy as np
from functools import lru_cache

# -*- coding: utf-8 -*-
# /**
#  * 各地图API坐标系统比较与转换;
#  * WGS84坐标系：即地球坐标系，国际上通用的坐标系。设备一般包含GPS芯片或者北斗芯片获取的经纬度为WGS84地理坐标系,
#  * 谷歌地图采用的是WGS84地理坐标系（中国范围除外）;
#  * GCJ02坐标系：即火星坐标系，是由中国国家测绘局制订的地理信息系统的坐标系统。由WGS84坐标系经加密后的坐标系。谷歌中国地图和搜搜中国地图采用的是GCJ02地理坐标系；
#  * 3BD09坐标系：即百度坐标系，GCJ02坐标系经加密后的坐标系;
#  */
import math
from decimal import *

class GLM(object):
    def __init__(self, fileName):
        self.fileName_ = fileName

    # 读取tif数据集
    def readTif(self):
        dataset = gdal.Open(self.fileName_)
        if dataset is None:
            print(self.fileName_+"文件无法打开")
        return dataset

    # 获取仿射矩阵信息
    def Getgeotrans(self):
        dataset = self.readTif()
        # dataset = GLM.readTif(self)

        im_width = dataset.RasterXSize  # 栅格矩阵的列数
        im_height = dataset.RasterYSize  # 栅格矩阵的行数
        im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵
        print("放射矩阵:",im_geotrans)
        im_proj = dataset.GetProjection()  # 地图投影信息
        im_data = dataset.ReadAsArray(0, 0, im_width, im_height)  # 将数据写成数组，对应栅格矩阵

        del dataset
        return im_proj, im_geotrans, im_data



    #像素坐标和地理坐标仿射变换
    def CoordTransf(self, Xpixel, Ypixel, GeoTransform): #wgs坐标
        XGeo = GeoTransform[0]+GeoTransform[1]*Xpixel+Ypixel*GeoTransform[2]
        YGeo = GeoTransform[3]+GeoTransform[4]*Xpixel+Ypixel*GeoTransform[5]
        return XGeo, YGeo

    # 写文件，写成tif
    def write_img(self, outpath, im_proj, im_geotrans, im_data):
        # gdal数据类型包括
        # gdal.GDT_Byte,
        # gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
        # gdal.GDT_Float32, gdal.GDT_Float64

        # 判断栅格数据的数据类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32

        # 判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape

        # 创建文件
        driver = gdal.GetDriverByName("GTiff")  # 数据类型必须有，因为要计算需要多大内存空间
        dataset = driver.Create(outpath, im_width, im_height, im_bands, datatype)

        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影

        if im_bands == 1:
            dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
        else:
            for i in range(im_bands):
                dataset.GetRasterBand(i + 1).WriteArray(im_data[i])

        del dataset


if __name__ == "__main__":
    filelist = glob("data_tif/cut_example.tif") #tif路径
    # image_list = []
    for file in filelist:
        print('[INFO] {}'.format(file))
        glm = GLM(file)
        proj, geotrans, data = glm.Getgeotrans()
        print(data)
        depth,height, width= data.shape
        print("depth,height,width:",depth,height,width)
        for y in range(0, height):
            for x in range(0, width):
                #longitude 经度 latitude 维度
                lon,lat = glm.CoordTransf(y, x, geotrans)
                #print('lat={}, lon={} is on land:'.format(lon, lat))
                is_on_land=is_land(lat,lon)
                #print('lat={}, lon={} is on land: {}'.format(lon, lat, is_on_land))
                if is_on_land==False:
                    data[:, y, x] = 0
                else:
                    data[:,y,x] = 255
                #print(data[:, y, x])
        print(data)
        #保存路径
        glm.write_img("fenli_data/cut_example_out.tif", proj, geotrans,data[:, :, :])  # 写数据
        #glm.write_img(os.path.join(out_dir, os.path.basename(file).split('.')[0]+'.tif'), proj, geotrans, data[:, :, :])  # 写数据
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

结果：（1m分辨率，可以看到，误差很大，没法用）

500m分辨率 MODIS效果就好多了

经纬度1秒误差在30m左右，1度🟰3600秒，如果遥感影像的分辨率大 比如（1m，或者2m）想实现海陆分离不能根据经纬度坐标了，这时候建议裁剪所需区域手动标记制作掩膜。

3 tif 转 jpg

import os,sys
import cv2
import numpy as np
from skimage import io#使用IO库读取tif图片


def tif_jpg_transform(file_path_name, bgr_savepath_name):
    img = io.imread(file_path_name)#读取文件名
    img = img / img.max()#使其所有值不大于一
    img = img * 255 - 0.001  # 减去0.001防止变成负整型
    img = img.astype(np.uint8)#强制转换成8位整型
    # img = np.array([img,img,img])
    # img = img.transpose(1,2,0)
    print(img.shape)  # 显示图片大小和深度
    b = img[:, :, 0]  # 读取蓝通道
    g = img[:, :, 1]  # 读取绿通道
    r = img[:, :, 2]  # 读取红通道
    bgr = cv2.merge([r, g, b])  # 通道拼接
    cv2.imwrite(bgr_savepath_name, bgr)#图片存储


tif_file_path = r'yourpath'  # 为tif图片的文件夹路径
tif_fileList = os.listdir(tif_file_path)
for tif_file in tif_fileList:
    file_path_name = tif_file_path + '/' + tif_file
    print(file_path_name)
    jpg_path = r'yourpath' + '/' + tif_file.split('.')[0] + '.jpg' #.jpg图片的保存路径
    print(jpg_path)
    tif_jpg_transform(file_path_name, jpg_path)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

4 jpg图像裁剪为指定分辨率

import cv2
import numpy as np
import os
pic_path = '/Users/xxx/PycharmProjects/py38_tf/遥感大作业/cut/222_out.jpg' # 分割的图片的位置
pic_target = '/Users/xxx/PycharmProjects/py38_tf/遥感大作业/cut/1024x1024/' # 分割后的图片保存的文件夹
if not os.path.exists(pic_target):  #判断是否存在文件夹如果不存在则创建为文件夹
    os.makedirs(pic_target)
#要分割后的尺寸
cut_width = 1024
cut_length = 1024
# 读取要分割的图片，以及其尺寸等数据
picture = cv2.imread(pic_path)
(width, length, depth) = picture.shape
print(width,length,depth)
# 预处理生成0矩阵
pic = np.zeros((cut_width, cut_length, depth))
# 计算可以划分的横纵的个数
num_width = int(width / cut_width)
num_length = int(length / cut_length)
# for循环迭代生成
for i in range(0, num_width):
    for j in range(0, num_length):
        pic = picture[i*cut_width : (i+1)*cut_width, j*cut_length : (j+1)*cut_length, :]
        result_path = pic_target + '{}_{}.jpg'.format(i+1, j+1)
        cv2.imwrite(result_path, pic)

print("done!!!")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

5 常见的几何变换扩充数据集（opencv实现）

# 程序功能：python+opencv实现数据增强，并保存图片
import numpy as np
import cv2
import os
#数据增强后的图片保存路劲，不要有中文路径！！！
org_path_images="D:/pycharm/MyCode/demo/data/training/images/"
save_path_images="D:/pycharm/MyCode/demo/data/training/DA_images/"
org_path_labels="D:/pycharm/MyCode/demo/data/training/labels/"
save_path_labels="D:/pycharm/MyCode/demo/data/training/DA_labels/"
def SDA(org_path,save_path):
    i=100
    print(os.listdir(org_path))
    for info in os.listdir(org_path):
        info1 = os.path.join(org_path, info)  # 将路径与文件名结合起来就是每个文件的完整路径
        print(info1)
        img = cv2.imread(info1)
        #print("img:",img)
        # if(img==None):
        #     print("None!!!")
        # 水平镜像
        i += 1
        h_flip = cv2.flip(img, 1)
        cv2.imwrite(save_path+str(i)+'.png', h_flip)
        # 垂直镜像
        i += 1
        v_flip = cv2.flip(img, 0)
        cv2.imwrite(save_path+str(i)+'.png', v_flip)
        # 水平垂直镜像
        i += 1
        hv_flip = cv2.flip(img, -1)
        #cv2.imshow("Flipped Horizontally & Vertically", hv_flip)
        cv2.imwrite(save_path+str(i)+'.png', hv_flip)
        # 90度旋转
        i += 1
        rows, cols = img.shape[:2]
        M = cv2.getRotationMatrix2D((cols / 2, rows / 2), 90, 1)
        dst_90 = cv2.warpAffine(img, M, (cols, rows))
        cv2.imwrite(save_path+str(i)+'.png', dst_90)
        # 70度旋转
        i += 1
        rows, cols = img.shape[:2]
        M = cv2.getRotationMatrix2D((cols / 2, rows / 2), 70, 1)
        dst_70 = cv2.warpAffine(img, M, (cols, rows))
        #cv2.imshow("dst_70", dst_70)
        cv2.imwrite(save_path+str(i)+'.png', dst_70)
        # 60度旋转
        i += 1
        rows, cols = img.shape[:2]
        M = cv2.getRotationMatrix2D((cols / 2, rows / 2), 60, 1)
        dst_60 = cv2.warpAffine(img, M, (cols, rows))
       # cv2.imshow("dst_60", dst_60)
        cv2.imwrite(save_path+str(i)+'.png', dst_60)
        # 50度旋转
        i += 1
        rows, cols = img.shape[:2]
        M = cv2.getRotationMatrix2D((cols / 2, rows / 2), 50, 1)
        dst_50 = cv2.warpAffine(img, M, (cols, rows))
        #cv2.imshow("dst_50", dst_50)
        cv2.imwrite(save_path+str(i)+'.png', dst_50)
        # 45度旋转
        i += 1
        rows, cols = img.shape[:2]
        M = cv2.getRotationMatrix2D((cols / 2, rows / 2), 45, 1)
        dst_45 = cv2.warpAffine(img, M, (cols, rows))
        #cv2.imshow("dst_45", dst_45)
        cv2.imwrite(save_path+str(i)+'.png', dst_45)
        # 40度旋转
        i += 1
        rows, cols = img.shape[:2]
        M = cv2.getRotationMatrix2D((cols / 2, rows / 2), 40, 1)
        dst_40 = cv2.warpAffine(img, M, (cols, rows))
        #cv2.imshow("dst_40", dst_40)
        cv2.imwrite(save_path+str(i)+'.png', dst_40)
        # 30度旋转
        i += 1
        rows, cols = img.shape[:2]
        M = cv2.getRotationMatrix2D((cols / 2, rows / 2), 30, 1)
        dst_30 = cv2.warpAffine(img, M, (cols, rows))
        #cv2.imshow("dst_30", dst_30)
        cv2.imwrite(save_path+str(i)+'.png', dst_30)
        # 20度旋转
        i += 1
        rows, cols = img.shape[:2]
        M = cv2.getRotationMatrix2D((cols / 2, rows / 2), 20, 1)
        dst_20 = cv2.warpAffine(img, M, (cols, rows))
        #cv2.imshow("dst_20", dst_20)
        cv2.imwrite(save_path+str(i)+'.png', dst_20)
        # 缩放
        # height, width = img.shape[:2]
        # res = cv2.resize(img, (2 * width, 2 * height))
        # cv2.imshow("large", res)
        # cv2.imwrite(save_path+info+'res.png', res)
        # 仿射变换
        # 对图像进行变换（三点得到一个变换矩阵）
        # 我们知道三点确定一个平面，我们也可以通过确定三个点的关系来得到转换矩阵
        # 然后再通过warpAffine来进行变换
        i += 1
        point1 = np.float32([[50, 50], [300, 50], [50, 200]])
        point2 = np.float32([[10, 100], [300, 50], [100, 250]])
        M = cv2.getAffineTransform(point1, point2)
        dst1 = cv2.warpAffine(img, M, (cols, rows), borderValue=(255, 255, 255))
        #cv2.imshow("affine transformation", dst1)
        cv2.imwrite(save_path+str(i)+'.png', dst1)
    cv2.waitKey(0)

if __name__=="__main__":
    #图片跟标签做同样的增强操作
    #print("ghi8gweohboier")
    SDA(org_path_images,save_path_images)
    SDA(org_path_labels,save_path_labels)


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112