visualize.py

# Copyright (c) 2022 PaddlePaddle Authors. All Rights Reserved. 
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import warnings

warnings.filterwarnings("ignore", category=DeprecationWarning)

import os.path as osp
from pathlib import Path
from typing import List, Tuple, Union, Optional
import webbrowser
import numpy as np
from folium import folium, Map, LayerControl
from folium.raster_layers import TileLayer, ImageOverlay
from paddlers.transforms.functions import to_uint8

try:
    from osgeo import gdal, osr
except:
    import gdal
    import osr

CHINA_TILES = (
    "GeoQ China Community",
    "GeoQ China Street",
    "AMAP China",
    "TencentMap China",
    "BaiduMaps China", )


def map_display(mask_path: str,
                img_path: Optional[str]=None,
                band_list: Union[List[int], Tuple[int, ...], None]=None,
                save_path: Optional[str]=None,
                tiles: str="GeoQ China Community") -> folium.Map:
    """
    Show mask (and original image) on an online map.

    Args:
    mask_path (str): Path of predicted or ground-truth masks.
    img_path (str|None, optional): Path of the original image. Defaults to None.
    band_list (list[int]|tuple[int]|None, optional): 
        Bands to select from the original image for display (the band index starts from 1). 
        If None, use all bands. Defaults to None.
    save_path (str, optional): Path of the .html file to save the visualization results. 
        In Jupyter Notebook environments, 
        leave `save_path` as None to display the result immediately in the notebook. 
        Defaults to None.
    tiles (str): Map tileset to use. Chosen from the following list:
            - "GeoQ China Community", "GeoQ China Street" (from http://www.geoq.cn/)
            - "AMAP China" (from https://www.amap.com/)
            - "TencentMap China" (from https://map.qq.com/)
            - "BaiduMaps China" (from https://map.baidu.com/)

        Defaults to "GeoQ China Community".

        * All tilesets have been corrected through public algorithms from the Internet.
        * Please read the relevant terms of use carefully:
            - GeoQ [GISUNI]  (http://geoq.cn/useragreement.html)
            - AMap [AutoNavi]  (https://wap.amap.com/doc/serviceitem.html)
            - Tencent Map  (https://ugc.map.qq.com/AppBox/Landlord/serveagreement.html)
            - Baidu Map  (https://map.baidu.com/zt/client/service/index.html)

    Returns:
        folium.Map: An example of folium map.
    """

    if tiles not in CHINA_TILES:
        raise ValueError("The `tiles` must in {}, not {}.".format(CHINA_TILES,
                                                                  tiles))
    fmap = Map(
        tiles=tiles,
        min_zoom=1,
        max_zoom=24, )
    if img_path is not None:
        layer, _ = Raster(img_path, band_list).get_layer()
        layer.add_to(fmap)
    layer, center = Raster(mask_path).get_layer()
    layer.add_to(fmap)
    if center is not None:
        fmap.location = center
    fmap.fit_bounds(layer.bounds)
    LayerControl().add_to(fmap)
    if save_path:
        fmap.save(save_path)
        webbrowser.open(save_path)
    return fmap


class OpenAsEPSG4326Error(Exception):
    pass


class Raster:
    def __init__(
            self,
            path: str,
            band_list: Union[List[int], Tuple[int, ...], None]=None) -> None:
        self.src_data = Converter.open_as_WGS84(path)
        if self.src_data is None:
            raise OpenAsEPSG4326Error("Faild to open {} in EPSG:4326.".format(
                path))
        self.name = Path(path).stem
        self.set_bands(band_list)
        self._get_info()

    def get_layer(self) -> ImageOverlay:
        layer = ImageOverlay(self._get_array(), self.wgs_range, name=self.name)
        return layer, self.wgs_center

    def set_bands(self,
                  band_list: Union[List[int], Tuple[int, ...], None]) -> None:
        self.bands = self.src_data.RasterCount
        if band_list is None:
            if self.bands == 3:
                band_list = [1, 2, 3]
            else:
                band_list = [1]
        band_list_lens = len(band_list)
        if band_list_lens not in (1, 3):
            raise ValueError("The lenght of band_list must be 1 or 3, not {}.".
                             format(str(band_list_lens)))
        if max(band_list) > self.bands or min(band_list) < 1:
            raise ValueError("The range of band_list must within [1, {}].".
                             format(str(self.bands)))
        self.band_list = band_list

    def _get_info(self) -> None:
        self.width = self.src_data.RasterXSize
        self.height = self.src_data.RasterYSize
        self.geotf = self.src_data.GetGeoTransform()
        self.proj = self.src_data.GetProjection()  # WGS84
        self.wgs_range = self._get_WGS84_range()
        self.wgs_center = self._get_WGS84_center()

    def _get_WGS84_range(self) -> List[List[float]]:
        converter = Converter(self.proj, self.geotf)
        lat1, lon1 = converter.xy2latlon(self.height - 1, 0)
        lat2, lon2 = converter.xy2latlon(0, self.width - 1)
        return [[lon1, lat1], [lon2, lat2]]

    def _get_WGS84_center(self) -> List[float]:
        clat = (self.wgs_range[0][0] + self.wgs_range[1][0]) / 2
        clon = (self.wgs_range[0][1] + self.wgs_range[1][1]) / 2
        return [clat, clon]

    def _get_array(self) -> np.ndarray:
        band_array = []
        for b in self.band_list:
            band_i = self.src_data.GetRasterBand(b).ReadAsArray()
            band_array.append(band_i)
        ima = np.stack(band_array, axis=0)
        if self.bands == 1:
            # the type is complex means this is a SAR data
            if isinstance(type(ima[0, 0]), complex):
                ima = abs(ima)
            ima = ima.squeeze()
        else:
            ima = ima.transpose((1, 2, 0))
        ima = to_uint8(ima, True)
        return ima


class Converter:
    def __init__(self, proj: str, geotf: tuple) -> None:
        # source data
        self.source = osr.SpatialReference()
        self.source.ImportFromWkt(proj)
        self.geotf = geotf
        # target data
        self.target = osr.SpatialReference()
        self.target.ImportFromEPSG(4326)

    @classmethod
    def open_as_WGS84(self, path: str) -> gdal.Dataset:
        if not osp.exists(path):
            raise FileNotFoundError("{} not found.".format(path))
        result = gdal.Warp("", path, dstSRS="EPSG:4326", format="VRT")
        return result

    def xy2latlon(self, row: int, col: int) -> List[float]:
        px = self.geotf[0] + col * self.geotf[1] + row * self.geotf[2]
        py = self.geotf[3] + col * self.geotf[4] + row * self.geotf[5]
        ct = osr.CoordinateTransformation(self.source, self.target)
        coords = ct.TransformPoint(px, py)
        return coords[:2]