image_tile.py 13.4 KB
# coding=utf-8
#author:        4N
#createtime:    2021/3/24
#email:         nheweijun@sina.com

from app.util import *
import traceback
from osgeo import gdal
from osgeo.gdal import *
from numpy import ndarray
import numpy
from flask import Response
import io
import os
from PIL import Image
from app.util.slice_scheme import slice_scheme
import time
import cv2
from .image_tile_center import create_by_opencv
def api(level,row,col):
    result = {}
    parameter: dict = get_parameter()


    try:

        bands = [1, 2, 3]

        image_list=[
            {"origin_extent":[111.604613312, 29.0171545762, 111.653989358, 29.0531633509],
             "path" : os.path.join(os.path.dirname(__file__), "data", "江南_03.tif"),
             "xysize":[47646,28201],
             "max_level":7},

            {"origin_extent": [111.705644552, 28.9864085959, 111.737115887, 29.079435995],
             "path": os.path.join(os.path.dirname(__file__), "data", "江南_01.tif"),
             "xysize":[30116, 73219],
             "max_level":8},

            {"origin_extent": [111.639350712, 28.9170588759, 111.751508603, 29.032941696],
             "path": os.path.join(os.path.dirname(__file__), "data", "江南_02.tif"),
             "xysize": [108488, 90777],
             "max_level": 9}
        ]



        slice_para = {'rows': 256.0, 'cols': 256.0, 'x': -180.0, 'y': 90.0, 'dpi': 96.0,
                      '0': {'scale': 590995186.11750006, 'resolution': 1.4062500000000004},
                      '1': {'scale': 295497593.05875003, 'resolution': 0.7031250000000002},
                      '2': {'scale': 147748796.52937502, 'resolution': 0.3515625000000001},
                      '3': {'scale': 73874398.26468751, 'resolution': 0.17578125000000006},
                      '4': {'scale': 36937199.132343754, 'resolution': 0.08789062500000003},
                      '5': {'scale': 18468599.566171877, 'resolution': 0.043945312500000014},
                      '6': {'scale': 9234299.783085939, 'resolution': 0.021972656250000007},
                      '7': {'scale': 4617149.891542969, 'resolution': 0.010986328125000003},
                      '8': {'scale': 2308574.9457714846, 'resolution': 0.005493164062500002},
                      '9': {'scale': 1154287.4728857423, 'resolution': 0.002746582031250001},
                      '10': {'scale': 577143.7364428712, 'resolution': 0.0013732910156250004},
                      '11': {'scale': 288571.8682214356, 'resolution': 0.0006866455078125002},
                      '12': {'scale': 144285.9341107178, 'resolution': 0.0003433227539062501},
                      '13': {'scale': 72142.9670553589, 'resolution': 0.00017166137695312505},
                      '14': {'scale': 36071.48352767945, 'resolution': 8.583068847656253e-05},
                      '15': {'scale': 18035.741763839724, 'resolution': 4.2915344238281264e-05},
                      '16': {'scale': 9017.870881919862, 'resolution': 2.1457672119140632e-05},
                      '17': {'scale': 4508.935440959931, 'resolution': 1.0728836059570316e-05},
                      '18': {'scale': 2254.4677204799655, 'resolution': 5.364418029785158e-06},
                      '19': {'scale': 1127.2338602399827, 'resolution': 2.682209014892579e-06},
                      '20': {'scale': 563.6169301199914, 'resolution': 1.3411045074462895e-06}}

        if parameter.get("leaflet"):
            slice_para = {'rows': 256.0, 'cols': 256.0, 'x': -180.0, 'y': 90.0, 'dpi': 96.0,
                          '0': {'scale': 295497593.05875003, 'resolution': 0.7031250000000002},
                          '1': {'scale': 147748796.52937502, 'resolution': 0.3515625000000001},
                          '2': {'scale': 73874398.26468751, 'resolution': 0.17578125000000006},
                          '3': {'scale': 36937199.132343754, 'resolution': 0.08789062500000003},
                          '4': {'scale': 18468599.566171877, 'resolution': 0.043945312500000014},
                          '5': {'scale': 9234299.783085939, 'resolution': 0.021972656250000007},
                          '6': {'scale': 4617149.891542969, 'resolution': 0.010986328125000003},
                          '7': {'scale': 2308574.9457714846, 'resolution': 0.005493164062500002},
                          '8': {'scale': 1154287.4728857423, 'resolution': 0.002746582031250001},
                          '9': {'scale': 577143.7364428712, 'resolution': 0.0013732910156250004},
                          '10': {'scale': 288571.8682214356, 'resolution': 0.0006866455078125002},
                          '11': {'scale': 144285.9341107178, 'resolution': 0.0003433227539062501},
                          '12': {'scale': 72142.9670553589, 'resolution': 0.00017166137695312505},
                          '13': {'scale': 36071.48352767945, 'resolution': 8.583068847656253e-05},
                          '14': {'scale': 18035.741763839724, 'resolution': 4.2915344238281264e-05},
                          '15': {'scale': 9017.870881919862, 'resolution': 2.1457672119140632e-05},
                          '16': {'scale': 4508.935440959931, 'resolution': 1.0728836059570316e-05},
                          '17': {'scale': 2254.4677204799655, 'resolution': 5.364418029785158e-06},
                          '18': {'scale': 1127.2338602399827, 'resolution': 2.682209014892579e-06},
                          '19': {'scale': 563.6169301199914, 'resolution': 1.3411045074462895e-06}}

        # 转换参数
        new_para ={}
        for key in parameter.keys():
            new_para[key.lower()] = parameter[key]
        parameter=new_para
        if parameter.get("tilematrix"):
            if parameter.get("tilematrix").__contains__(":"):
                level = int(parameter.get("tilematrix").split(":")[-1])
            else:
                level = int(parameter.get("tilematrix"))
        if parameter.get("tilerow"):
            row = int(parameter.get("tilerow"))
        if parameter.get("tilecol"):
            col = int(parameter.get("tilecol"))

        image_type = parameter.get("format") if parameter.get("format") else "image/png"

        quality = int(parameter.get("quality")) if parameter.get("quality") else 30

        extent = slice_scheme.get_polygon(slice_para, level, row, col)


        pixel_array = numpy.zeros((256, 256,3), dtype=int)
        ceng = 2
        for band in bands:

            empty = numpy.zeros((256, 256), dtype=int)+65536

            for im in image_list:

                # 自决定金字塔等级
                xysize = im.get("xysize")
                origin_extent = im.get("origin_extent")
                max_level = im.get("max_level")

                # 超出空间范围
                if extent[2]<origin_extent[0] or extent[0]>origin_extent[2] or extent[1]>origin_extent[3] or extent[3]<origin_extent[1]:
                    pass
                # 空间范围相交
                else:
                    image_level = determine_level(xysize, origin_extent, extent, max_level)
                    path = im.get("path")

                    # print(image_level)

                    image: Dataset = gdal.Open(path, 0)

                    band_data: Band = image.GetRasterBand(band)

                    if image_level == -1:
                        overview = band_data
                    else:
                        try:
                            overview: Band = band_data.GetOverview(image_level)
                        except:
                            raise Exception("该影像不存在该级别的金字塔数据!")
                    ox = overview.XSize
                    oy = overview.YSize

                    # 网格大小
                    grid_x = (origin_extent[2] - origin_extent[0]) / (ox * 1.0)
                    grid_y = (origin_extent[3] - origin_extent[1]) / (oy * 1.0)

                    # 完全在影像范围内
                    if extent[0]>origin_extent[0] and extent[1]>origin_extent[1] and extent[2]<origin_extent[2] and extent[3]<origin_extent[3]:

                        t1 = time.time()
                        # 网格偏移量

                        off_x = math.floor((extent[0]-origin_extent[0])/grid_x)
                        off_y = math.floor((origin_extent[3] -extent[3]) / grid_y)

                        # 截取后网格个数
                        x_g = math.ceil((extent[2]-extent[0])/grid_x)

                        y_g= math.ceil((extent[3]-extent[1])/grid_y)
                        t2 = time.time()
                        # print(t2-t1)
                        overview_raster:ndarray = overview.ReadAsArray(off_x,off_y,x_g,y_g,256,256)
                        t3 = time.time()
                        # print(t3-t2)

                        mask1 = numpy.zeros((256, 256), dtype=int)
                        mask2 = numpy.zeros((256, 256), dtype=int)
                        mask1[overview_raster == 65536] = 1
                        mask2[overview_raster != 65536] = 1

                        empty = empty*mask1+overview_raster*mask2
                        t4 = time.time()
                        # print(t4-t3)
                    # 部分相交
                    else:

                        inter_extent = [0,0,0,0]
                        inter_extent[0] = origin_extent[0] if origin_extent[0]>extent[0] else extent[0]
                        inter_extent[1] = origin_extent[1] if origin_extent[1] > extent[1] else extent[1]
                        inter_extent[2] = origin_extent[2] if origin_extent[2] < extent[2] else extent[2]
                        inter_extent[3] = origin_extent[3] if origin_extent[3] < extent[3] else extent[3]

                        # 网格偏移量
                        off_x = math.floor((inter_extent[0]-origin_extent[0])/grid_x)
                        off_y = math.floor((origin_extent[3] -inter_extent[3]) / grid_y)

                        # 截取后网格个数
                        x_g = math.floor((inter_extent[2]-inter_extent[0])/grid_x)
                        y_g= math.floor((inter_extent[3]-inter_extent[1])/grid_y)

                        # 相对于出图的偏移量

                        #出图的网格大小
                        out_grid_x = (extent[2] - extent[0]) / (256 * 1.0)
                        out_grid_y = (extent[3] - extent[1]) / (256 * 1.0)

                        out_off_x = int(math.ceil((inter_extent[0]-extent[0])/out_grid_x))
                        out_off_y = int(math.ceil((extent[3] -inter_extent[3]) / out_grid_y))

                        out_x_g = int(math.floor((inter_extent[2]-inter_extent[0])/out_grid_x))
                        out_y_g= int(math.floor((inter_extent[3]-inter_extent[1])/out_grid_y))

                        # 相交部分在出图的哪个位置

                        overview_raster:ndarray = overview.ReadAsArray(off_x,off_y,x_g,y_g,out_x_g,out_y_g)

                        mask1 = numpy.zeros((out_y_g,out_x_g), dtype=int)

                        mask2 = numpy.zeros((out_y_g,out_x_g), dtype=int)

                        mask1[overview_raster == 65536] = 1
                        mask2[overview_raster != 65536] = 1

                        empty[out_off_y:out_off_y + out_y_g, out_off_x:out_off_x + out_x_g] = empty[out_off_y:out_off_y + out_y_g,out_off_x:out_off_x + out_x_g]*mask1 + overview_raster*mask2

                    # 关闭句柄
                    del image



            # opencv 的颜色排列为GBR
            pixel_array[:,:,ceng]=empty
            ceng-=1

        t5 = time.time()
        # print(t4-t3)

        #将图片生成在内存中,然后直接返回response
        im_data = create_by_opencv(image_type, pixel_array, quality)

        return Response(im_data, mimetype=image_type.lower())

    except Exception as e:
        print(traceback.format_exc())
        result["state"] = -1
        result["message"] = e.__str__()
        return result





def determine_level(xysize,origin_extent,extent,max_level):
    x = xysize[0]
    y = xysize[1]

    level = -1
    pixel = x*y * (((extent[2]-extent[0])*(extent[3]-extent[1]))/((origin_extent[2]-origin_extent[0])*(origin_extent[3]-origin_extent[1])))

    while pixel>100000 and level<max_level-1:
        level+=1
        x=x/2
        y=y/2
        pixel = x * y * (((extent[2] - extent[0]) * (extent[3] - extent[1])) / (
                (origin_extent[2] - origin_extent[0]) * (origin_extent[3] - origin_extent[1])))
    return level



def create_by_pil(image_type,pixel_list,quality):
    buffer = io.BytesIO()
    if image_type.__eq__("image/jpeg") or image_type.__eq__("image/jpg"):
        im_type = "jpeg"
        data = list(zip(pixel_list[0][0], pixel_list[1][0], pixel_list[2][0]))
        image_out = Image.new("RGB", (256, 256))

    else:
        im_type = "png"
        four = [0 if x.__eq__(65536) else 255 for x in pixel_list[0][0]]
        data = list(zip(pixel_list[0][0], pixel_list[1][0], pixel_list[2][0], four))
        image_out = Image.new("RGBA", (256, 256))
    t6 = time.time()

    image_out.putdata(data)
    image_out.save(buffer, im_type, quality=quality, optimize=True)
    im_data = buffer.getvalue()
    buffer.close()
    t7 = time.time()
    return im_data


api_doc={
"tags":["影像接口"],
"parameters":[

],
"responses":{
    200:{
        "schema":{
            "properties":{
            }
        }
        }
    }
}