"""pyproj库未使用"""
import pyproj
import shapely
from shapely import wkt
import math

"""wgs84_to_webmercator函数未使用"""
def wgs84_to_webmercator(lat_in,lon_in):
  input_crs = pyproj.CRS.from_string("espg:4326")
  output_crs = pyproj.CRS.from_string("espg:3857")
  coord_trans = pyproj.Transformer.from_crs(input_crs,output_crs)
  lat_out,lon_out = coord_trans.transform(lat_in,lon_in)
  return lat_out,lon_out

def deg2num(lat_deg,lon_deg,zoom_level):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom_level
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.asinh(math.tan(lat_deg)) / math.pi) / 2.0 * n)
  return xtile,ytile

def num2deg(xtile,ytile,zoom_level):
  n = 2.0 ** zoom_level
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)

  lon_deg_plus = (xtile + 1) / n * 360.0 - 180.0
  lat_rad_plus = math.atan(math.sinh(math.pi * (1 - 2 * (ytile + 1) / n)))
  lat_deg_plus = math.degrees(lat_rad_plus)

  return lat_deg,lon_deg,lat_deg_plus,lon_deg_plus

def get_tile_list(area_geom_wkt,zoom_level):
  tile_list = []
  area_geometry = wkt.loads(area_geom_wkt)
  """获取范围框的矩形外包围框坐标"""
  (x_min,y_min,x_max,y_max) = area_geometry.bounds

  """获取矩形外包围框坐标涉及到瓦片编号的x、y的最大最小值"""
  start_x,end_y = deg2num(y_min,x_min,zoom_level)
  end_x,start_y = deg2num(y_max,x_max,zoom_level)

  """遍历矩形外包围框坐标涉及到的全部瓦片编号，判断与原始范围框是否存在相交关系，是则正确"""
  for x_tile in range(start_x,end_x + 1):
     for y_tile in range(start_y,end_y + 1):
        """根据瓦片编号，获取瓦片对应的角点坐标"""
        lat_deg,lon_deg,lat_deg_plus,lon_deg_plus = num2deg(x_tile,y_tile,zoom_level)
        """将瓦片角点坐标转换为矩形几何"""
        bbox = shapely.geometry.box(lon_deg_plus,lat_deg,lon_deg,lat_deg_plus)
        """判断瓦片范围与原始范围框的相交关系"""
        if area_geometry.intersects(bbox):
          """去重处理"""
          if '_'.join([str(x_tile),str(y_tile)]) not in tile_list:
             tile_list.append('_'.join([str(x_tile),str(y_tile)]))
          else:
             pass
  return tile_list