# 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]