dms-client/scan_data/example.py

from xml.dom import minidom
from osgeo import gdal
from PIL import Image
import numpy as np
import os


def uint16to8(bands, lower_percent=0.001, higher_percent=99.999):
    """
    拉伸图像：图片16位转8位
    :param bands: 输入栅格数据
    :param lower_percent: 最低百分比
    :param higher_percent: 最高百分比
    :return:
    """
    out = np.zeros_like(bands, dtype=np.uint8)
    n = bands.shape[0]
    for i in range(n):
        a = 0  # np.min(band)
        b = 255  # np.max(band)
        c = np.percentile(bands[i, :, :], lower_percent)
        d = np.percentile(bands[i, :, :], higher_percent)
        t = a + (bands[i, :, :] - c) * (b - a) / (d - c)
        t[t < a] = a
        t[t > b] = b
        out[i, :, :] = t
    return out


def createXML(metadata, xlm_file):
    """
    创建xlm文件并写入字典
    :param metadata: 元数据信息
    :param xlm_file: xlm文件
    :return:
    """
    # 创建一个空的文档
    document = minidom.Document()  # 创建DOM文档对象
    # 创建一个根节点对象
    root = document.createElement('ProductMetaData')
    # 设置根节点的属性
    # root.setAttribute('', '')
    # 将根节点添加到文档对象中
    document.appendChild(root)
    # 字典转xml
    for key in metadata:
        # 创建父节点
        node_name = document.createElement(key)
        # 给父节点设置文本
        node_name.appendChild(document.createTextNode(str(metadata[key])))
        # 将各父节点添加到根节点
        root.appendChild(node_name)
    # 写入xlm文档
    with open(xlm_file, 'w', encoding='utf-8') as f:
        document.writexml(f, indent='\t', newl='\n', addindent='\t', encoding='utf-8')
    f.close()


def GetJPSSData(in_file, xml_path, thumbnail_path):
    """
    获取联合极轨卫星系统（JPSS-1）元数据：NOAA-20(Joint Polar Satellite System spacecraft)
    :param xml_path:
    :param thumbnail_path:
    :param in_file:
    :return: 元数据字典
    """
    try:
        # 生成缩略图
        in_path, basename = os.path.split(in_file)
        ThumbnailName = os.path.splitext(basename)[0] + "_thumb.jpg"
        ThumbnailPath = os.path.join(thumbnail_path, ThumbnailName)

        gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
        in_datasets = gdal.Open(in_file)
        meta_data = in_datasets.GetMetadata()
        # 取出子数据集
        datasets = in_datasets.GetSubDatasets()
        red_data = gdal.Open(datasets[0][0]).ReadAsArray()
        nir_data = gdal.Open(datasets[3][0]).ReadAsArray()
        swir_data = gdal.Open(datasets[9][0]).ReadAsArray()
        img_data = np.array([red_data, nir_data, swir_data])
        img_data = uint16to8(img_data)
        # Array转Image
        img_data2 = np.transpose(img_data, (1, 2, 0))
        img_data2 = img_data2[:, :, ::-1]
        img = Image.fromarray(img_data2)
        # 压缩图片大小
        if img_data.shape[1] > img_data.shape[2]:
            width = 512
            height = int(width / img_data.shape[1] * img_data.shape[2])
        else:
            height = 512
            width = int(height / img_data.shape[1] * img_data.shape[2])
        img.thumbnail((width, height))
        img.save(ThumbnailPath, "PNG")

        # 释放内存
        del in_datasets
        del img_data
        del img_data2
        del img

        # 生成XML文件
        xmlFileName = os.path.splitext(basename)[0] + ".xml"
        xmlPath = os.path.join(xml_path, xmlFileName)
        createXML(meta_data, xmlPath)

        # 产品日期
        ProductionTime = meta_data['ProductionTime']
        StartTime = meta_data['StartTime']
        EndTime = meta_data['EndTime']

        #  其他信息
        ImageGSD = str(meta_data['LongName']).split(" ")[-1]
        Bands = str(meta_data['title']).split(" ")[1]

        # 中心经纬度
        productUpperLeftLat = meta_data['NorthBoundingCoordinate']  # 左上纬度
        productUpperLeftLong = meta_data['WestBoundingCoordinate']  # 左上经度
        productUpperRightLat = meta_data['NorthBoundingCoordinate']  # 右上纬度
        productUpperRightLong = meta_data['EastBoundingCoordinate']  # 右上经度
        productLowerLeftLat = meta_data['SouthBoundingCoordinate']  # 左下纬度
        productLowerLeftLong = meta_data['WestBoundingCoordinate']  # 左下经度
        productLowerRightLat = meta_data['SouthBoundingCoordinate']  # 右下纬度
        productLowerRightLong = meta_data['EastBoundingCoordinate']  # 右下纬度

        # 边界几何
        boundaryGeomStr = f'POLYGON(({productUpperLeftLong} {productUpperLeftLat},' \
                          f'{productUpperRightLong} {productUpperRightLat},' \
                          f'{productLowerRightLong} {productLowerRightLat},' \
                          f'{productLowerLeftLong} {productLowerLeftLat},' \
                          f'{productUpperLeftLong} {productUpperLeftLat}))'

        # 构建字典
        jpss_dict = {"ProduceTime": ProductionTime,
                     "StartTime": StartTime,
                     "EndTime": EndTime,
                     "CloudPercent": "",
                     # "TopLeftLatitude": productUpperLeftLat,
                     # "TopLeftLongitude": productUpperLeftLong,
                     # "TopRightLatitude": productUpperRightLat,
                     # "TopRightLongitude": productUpperRightLong,
                     # "BottomLeftLatitude": productLowerLeftLat,
                     # "BottomLeftLongitude": productLowerLeftLong,
                     # "BottomRightLatitude": productLowerRightLat,
                     # "BottomRightLongitude": productLowerRightLong,
                     "boundaryGeomStr": boundaryGeomStr,
                     "bands": Bands,
                     "ImageGSD": ImageGSD,
                     "ProjectedCoordinates": "",
                     "CollectionCode": "",
                     "ThumbnailPath": ThumbnailPath,
                     "ThumbnailName": ThumbnailName,
                     "xmlPath": xmlPath,
                     "xmlFileName": xmlFileName,
                     "DirectoryDepth": "day"}

        # 判断字典是否为空
        if not jpss_dict:
            return {"code": -1, "msg": "没有满足条件的数据字典..."}
        print(jpss_dict)
        return jpss_dict
    except Exception as e:
        print(str(e))
        return {"code": -1, "msg": str(e)}