GB/T 13989-2012国家基本比例尺地形图分幅和编号的实现

国家基本比例尺地形图分幅和编号的实现

1. js实现

/**
 1. 采用标准：国家基本比例尺地形图分幅和编号《GB/T 13989-2012》
 2. 1、根据经纬度计算图幅号
 3. 2、根据图幅号图幅范围计算
 4. 3、定义图幅信息
 */
export default class MapFrameUtils {
  /**
 5. 根据经纬度计算100万比例尺的行列号
 6. @param lon
 7. @param lat
 8. @return {{col: number, row: number, code: string}}
   */
  static calcMapFrameNumber100 = (lon, lat) => {
    let a = Math.trunc(lat / 4.0) + 1;
    let b = Math.trunc(lon / 6.0) + 31;
    let mateInfo = MapFrameUtils.mateData100.find(
      (item) => item.rowNumber === a
    );
    return {
      row: a,
      col: b,
      code: mateInfo.rowCode,
    };
  };
  /**
 9. 根据比例尺，经纬度计算地图的行列号
 10. @param scale
 11. @param lon
 12. @param lat
 13. @return {{col: number, code: *, row: number, digit: (number|*)}|null}
   */
  static calcMapFrameNumberNot100 = (scale, lon, lat) => {
    let mateData = MapFrameUtils.mateData.find((item) => item.scale === scale);
    if (mateData) {
      let cStep = ((lat % 4) / mateData.latD).toPrecision(4);
      let c = Math.round(4 / mateData.latD - Math.trunc(Number(cStep)));
      let dStep = ((lon % 6) / mateData.lonD).toPrecision(4);
      let d = Math.trunc(Number(dStep)) + 1;
      return {
        row: c,
        col: d,
        code: mateData.code,
        digit: mateData.digit,
      };
    } else {
      return null;
    }
  };
  /**
 14. 根据计算出来的行列号，比例尺，计算该格网在地图商的范围
 15. @param frame100
 16. @param frameNot100
 17. @param scale
 18. @return {{minY: number, minX: number, maxY: number, maxX: number}}
   */
  static calcMapFrameExtent = (frame100, frameNot100, scale) => {
    let mateData = MapFrameUtils.mateData.find((item) => item.scale === scale);
    let lon = (frame100.col - 31) * 6 + (frameNot100.col - 1) * mateData.lonD;
    let lat =
      (frame100.row - 1) * 4 +
      (4 / mateData.latD - frameNot100.row) * mateData.latD;
 
    return {
      minX: lon,
      minY: lat,
      maxX: lon + mateData.lonD,
      maxY: lat + mateData.latD,
    };
  };
  /**
 19. 根据比例尺和经纬度，获取所在图幅范围的编号、行列号以及范围
 20. @param scale
 21. @param lon
 22. @param lat
 23. @return {null|{extent: {minY: number, minX: number, maxY: number, maxX: number}, col: number, code: string, row: number}}
   */
  static calcMapFrameNumber = (scale, lon, lat) => {
    let mapFrameNumber100 = MapFrameUtils.calcMapFrameNumber100(lon, lat);
    let mapFrameNumberNot100 = MapFrameUtils.calcMapFrameNumberNot100(
      scale,
      lon,
      lat
    );
    if (mapFrameNumberNot100) {
      let code =
        mapFrameNumber100.code +
        "" +
        mapFrameNumber100.col +
        "" +
        mapFrameNumberNot100.code +
        "" +
        mapFrameNumberNot100.row
          .toString()
          .padStart(mapFrameNumberNot100.digit, "0") +
        "" +
        mapFrameNumberNot100.col
          .toString()
          .padStart(mapFrameNumberNot100.digit, "0");
      let extent = MapFrameUtils.calcMapFrameExtent(
        mapFrameNumber100,
        mapFrameNumberNot100,
        scale
      );
      return {
        code: code,
        extent: extent,
        row: mapFrameNumberNot100.row,
        col: mapFrameNumberNot100.col,
      };
    } else {
      return null;
    }
  };
  /**
 24. 100万比例尺下的行号对应的代码
   */
  static mateData100 = [
    {
      rowCode: "A",
      rowNumber: 1,
    },
    {
      rowCode: "B",
      rowNumber: 2,
    },
    {
      rowCode: "C",
      rowNumber: 3,
    },
    {
      rowCode: "D",
      rowNumber: 4,
    },
    {
      rowCode: "E",
      rowNumber: 5,
    },
    {
      rowCode: "F",
      rowNumber: 6,
    },
    {
      rowCode: "G",
      rowNumber: 7,
    },
    {
      rowCode: "H",
      rowNumber: 8,
    },
    {
      rowCode: "I",
      rowNumber: 9,
    },
    {
      rowCode: "J",
      rowNumber: 10,
    },
    {
      rowCode: "K",
      rowNumber: 11,
    },
    {
      rowCode: "L",
      rowNumber: 12,
    },
    {
      rowCode: "M",
      rowNumber: 13,
    },
    {
      rowCode: "N",
      rowNumber: 14,
    },
    {
      rowCode: "O",
      rowNumber: 15,
    },
  ];
  /**
 25. 非100万情况下，不同比例尺的原始数据
   */
  static mateData = [
    {
      scale: 500000,
      label: "1:50万",
      lonD: 3,
      latD: 2,
      code: "B",
      digit: 3,
    },
    {
      scale: 250000,
      label: "1:25万",
      lonD: 1.5,
      latD: 1,
      code: "C",
      digit: 3,
    },
    {
      scale: 100000,
      label: "1:10万",
      lonD: 0.5,
      latD: 0.33333333,
      code: "D",
      digit: 3,
    },
    {
      scale: 50000,
      label: "1:5万",
      lonD: 0.25,
      latD: 0.16666667,
      code: "E",
      digit: 3,
    },
    {
      scale: 25000,
      label: "1:2.5万",
      lonD: 0.125,
      latD: 0.08333333,
      code: "F",
      digit: 3,
    },
    {
      scale: 10000,
      label: "1:1万",
      lonD: 0.0625,
      latD: 0.04166667,
      code: "G",
      digit: 3,
    },
    {
      scale: 5000,
      label: "1:5000",
      lonD: 0.03125,
      latD: 0.02083333,
      code: "H",
      digit: 3,
    },
    {
      scale: 2000,
      label: "1:2000",
      lonD: 0.01041667,
      latD: 0.006944444,
      code: "I",
      digit: 3,
    },
    {
      scale: 1000,
      label: "1:1000",
      lonD: 0.005208333,
      latD: 0.003472222,
      code: "J",
      digit: 4,
    },
    {
      scale: 500,
      label: "1:500",
      lonD: 0.002603889,
      latD: 0.001736111,
      code: "K",
      digit: 4,
    },
  ];
}

26. python实现

import math

class MapFrameUtils:
    # 100万比例尺下的行号对应的代码
    mate_data_100 = [
        {"row_code": "A", "row_number": 1},
        {"row_code": "B", "row_number": 2},
        {"row_code": "C", "row_number": 3},
        {"row_code": "D", "row_number": 4},
        {"row_code": "E", "row_number": 5},
        {"row_code": "F", "row_number": 6},
        {"row_code": "G", "row_number": 7},
        {"row_code": "H", "row_number": 8},
        {"row_code": "I", "row_number": 9},
        {"row_code": "J", "row_number": 10},
        {"row_code": "K", "row_number": 11},
        {"row_code": "L", "row_number": 12},
        {"row_code": "M", "row_number": 13},
        {"row_code": "N", "row_number": 14},
        {"row_code": "O", "row_number": 15},
    ]

    # 非100万情况下，不同比例尺的原始数据
    mate_data = [
        {"scale": 500000, "label": "1:50万", "lon_d": 3, "lat_d": 2, "code": "B", "digit": 3},
        {"scale": 250000, "label": "1:25万", "lon_d": 1.5, "lat_d": 1, "code": "C", "digit": 3},
        {"scale": 100000, "label": "1:10万", "lon_d": 0.5, "lat_d": 0.33333333, "code": "D", "digit": 3},
        {"scale": 50000, "label": "1:5万", "lon_d": 0.25, "lat_d": 0.16666667, "code": "E", "digit": 3},
        {"scale": 25000, "label": "1:2.5万", "lon_d": 0.125, "lat_d": 0.08333333, "code": "F", "digit": 3},
        {"scale": 10000, "label": "1:1万", "lon_d": 0.0625, "lat_d": 0.04166667, "code": "G", "digit": 3},
        {"scale": 5000, "label": "1:5000", "lon_d": 0.03125, "lat_d": 0.02083333, "code": "H", "digit": 3},
        {"scale": 2000, "label": "1:2000", "lon_d": 0.01041667, "lat_d": 0.006944444, "code": "I", "digit": 3},
        {"scale": 1000, "label": "1:1000", "lon_d": 0.005208333, "lat_d": 0.003472222, "code": "J", "digit": 4},
        {"scale": 500, "label": "1:500", "lon_d": 0.002603889, "lat_d": 0.001736111, "code": "K", "digit": 4},
    ]

    @staticmethod
    def calc_map_frame_number_100(lon, lat):
        """
        根据经纬度计算100万比例尺的行列号
        :param lon: 经度
        :param lat: 纬度
        :return: dict containing row, col, and code
        """
        a = math.trunc(lat / 4.0) + 1
        b = math.trunc(lon / 6.0) + 31
        mate_info = next((item for item in MapFrameUtils.mate_data_100 if item["row_number"] == a), None)
        return {
            "row": a,
            "col": b,
            "code": mate_info["row_code"] if mate_info else ""
        }

    @staticmethod
    def calc_map_frame_number_not_100(scale, lon, lat):
        """
        根据比例尺，经纬度计算地图的行列号
        :param scale: 比例尺
        :param lon: 经度
        :param lat: 纬度
        :return: dict containing row, col, code, and digit, or None
        """
        mate_data = next((item for item in MapFrameUtils.mate_data if item["scale"] == scale), None)
        if mate_data:
            c_step = float(f"{(lat % 4) / mate_data['lat_d']:.4f}")
            c = round(4 / mate_data["lat_d"] - math.trunc(c_step))
            d_step = float(f"{(lon % 6) / mate_data['lon_d']:.4f}")
            d = math.trunc(d_step) + 1
            return {
                "row": c,
                "col": d,
                "code": mate_data["code"],
                "digit": mate_data["digit"]
            }
        return None

    @staticmethod
    def calc_map_frame_extent(frame_100, frame_not_100, scale):
        """
        根据计算出来的行列号，比例尺，计算该格网在地图上的范围
        :param frame_100: 100万比例尺的行列号
        :param frame_not_100: 非100万比例尺的行列号
        :param scale: 比例尺
        :return: dict containing minX, minY, maxX, maxY
        """
        mate_data = next((item for item in MapFrameUtils.mate_data if item["scale"] == scale), None)
        if mate_data:
            lon = (frame_100["col"] - 31) * 6 + (frame_not_100["col"] - 1) * mate_data["lon_d"]
            lat = (frame_100["row"] - 1) * 4 + (4 / mate_data["lat_d"] - frame_not_100["row"]) * mate_data["lat_d"]
            return {
                "min_x": lon,
                "min_y": lat,
                "max_x": lon + mate_data["lon_d"],
                "max_y": lat + mate_data["lat_d"]
            }
        return {}

    @staticmethod
    def calc_map_frame_number(scale, lon, lat):
        """
        根据比例尺和经纬度，获取所在图幅范围的编号、行列号以及范围
        :param scale: 比例尺
        :param lon: 经度
        :param lat: 纬度
        :return: dict containing code, extent, row, and col, or None
        """
        map_frame_number_100 = MapFrameUtils.calc_map_frame_number_100(lon, lat)
        map_frame_number_not_100 = MapFrameUtils.calc_map_frame_number_not_100(scale, lon, lat)
        
        if map_frame_number_not_100:
            code = (
                f"{map_frame_number_100['code']}"
                f"{map_frame_number_100['col']}"
                f"{map_frame_number_not_100['code']}"
                f"{str(map_frame_number_not_100['row']).zfill(map_frame_number_not_100['digit'])}"
                f"{str(map_frame_number_not_100['col']).zfill(map_frame_number_not_100['digit'])}"
            )
            extent = MapFrameUtils.calc_map_frame_extent(
                map_frame_number_100,
                map_frame_number_not_100,
                scale
            )
            return {
                "code": code,
                "extent": extent,
                "row": map_frame_number_not_100["row"],
                "col": map_frame_number_not_100["col"]
            }
        return None

    @staticmethod
    def calc_bbox_from_code(code):
        """
        根据图幅编号计算边界框（bbox）
        :param code: 图幅编号，例如 "A31B001001"
        :return: dict containing min_x, min_y, max_x, max_y, or None if invalid
        """
        if not code or not isinstance(code, str):
            return None

        # Step 1: Extract row code (first character)
        row_code = code[0]
        mate_info_100 = next((item for item in MapFrameUtils.mate_data_100 if item["row_code"] == row_code), None)
        if not mate_info_100:
            return None
        row_100 = mate_info_100["row_number"]

        # Step 2: Extract column number for 1:1,000,000 scale
        # The column number follows the row code until the scale-specific code
        i = 1
        col_str = ""
        while i < len(code) and code[i].isdigit():
            col_str += code[i]
            i += 1
        try:
            col_100 = int(col_str)
        except ValueError:
            return None

        # Step 3: Extract scale-specific code (single character after column number)
        if i >= len(code):
            return None
        scale_code = code[i]
        mate_data = next((item for item in MapFrameUtils.mate_data if item["code"] == scale_code), None)
        if not mate_data:
            return None
        scale = mate_data["scale"]
        digit = mate_data["digit"]

        # Step 4: Extract scale-specific row and column numbers
        remaining = code[i + 1:]
        if len(remaining) != 2 * digit:
            return None
        try:
            row_not_100 = int(remaining[:digit])
            col_not_100 = int(remaining[digit:])
        except ValueError:
            return None

        # Step 5: Construct frame_100 and frame_not_100
        frame_100 = {"row": row_100, "col": col_100}
        frame_not_100 = {"row": row_not_100, "col": col_not_100}

        # Step 6: Calculate extent (bbox)
        extent = MapFrameUtils.calc_map_frame_extent(frame_100, frame_not_100, scale)
        if not extent:
            return None

        return {
            "min_x": extent["min_x"],
            "min_y": extent["min_y"],
            "max_x": extent["max_x"],
            "max_y": extent["max_y"]
        }

    @staticmethod
    def calc_map_frames_in_range(min_x, min_y, max_x, max_y, scale):
        """
        根据给定的坐标范围和比例尺，计算所有图幅编号及边界框
        :param min_x: 最小经度
        :param min_y: 最小纬度
        :param max_x: 最大经度
        :param max_y: 最大纬度
        :param scale: 比例尺
        :return: list of dicts containing code and bbox (min_x, min_y, max_x, max_y)
        """
        result = []

        # Validate scale
        mate_data = next((item for item in MapFrameUtils.mate_data if item["scale"] == scale), None)
        if not mate_data:
            return result
        lon_d = mate_data["lon_d"]
        lat_d = mate_data["lat_d"]
        digit = mate_data["digit"]
        scale_code = mate_data["code"]

        # Step 1: Determine 1:1,000,000 scale grid range
        min_row_100 = math.trunc(min_y / 4.0) + 1
        max_row_100 = math.trunc(max_y / 4.0) + 1
        min_col_100 = math.trunc(min_x / 6.0) + 31
        max_col_100 = math.trunc(max_x / 6.0) + 31

        # Step 2: Iterate over 1:1,000,000 scale grid cells
        for row_100 in range(min_row_100, max_row_100 + 1):
            mate_info_100 = next((item for item in MapFrameUtils.mate_data_100 if item["row_number"] == row_100), None)
            if not mate_info_100:
                continue
            row_code = mate_info_100["row_code"]

            for col_100 in range(min_col_100, max_col_100 + 1):
                # Calculate 1:1,000,000 scale grid boundaries
                lon_start = (col_100 - 31) * 6
                lat_start = (row_100 - 1) * 4
                lon_end = lon_start + 6
                lat_end = lat_start + 4

                # Step 3: Determine scale-specific grid range within this 1:1,000,000 cell
                max_row_not_100 = int(4 / lat_d)  # Number of rows in scale-specific grid
                max_col_not_100 = int(6 / lon_d)  # Number of columns in scale-specific grid

                # Calculate scale-specific row and column ranges that intersect with input bbox
                min_row_not_100 = max(1, int((4 / lat_d) - math.ceil((lat_end - min_y) / lat_d)))
                max_row_not_100 = min(max_row_not_100, int((4 / lat_d) - math.floor((lat_end - max_y) / lat_d)))
                min_col_not_100 = max(1, math.floor((min_x - lon_start) / lon_d) + 1)
                max_col_not_100 = min(max_col_not_100, math.ceil((max_x - lon_start) / lon_d))

                # Step 4: Generate codes and bboxes for scale-specific grid cells
                for row_not_100 in range(min_row_not_100, max_row_not_100 + 1):
                    for col_not_100 in range(min_col_not_100, max_col_not_100 + 1):
                        # Construct frame_100 and frame_not_100
                        frame_100 = {"row": row_100, "col": col_100}
                        frame_not_100 = {"row": row_not_100, "col": col_not_100}

                        # Calculate extent
                        extent = MapFrameUtils.calc_map_frame_extent(frame_100, frame_not_100, scale)
                        if not extent:
                            continue

                        # Generate code
                        code = (
                            f"{row_code}"
                            f"{col_100}"
                            f"{scale_code}"
                            f"{str(row_not_100).zfill(digit)}"
                            f"{str(col_not_100).zfill(digit)}"
                        )

                        # Add to result
                        result.append({
                            "code": code,
                            "bbox": {
                                "min_x": extent["min_x"],
                                "min_y": extent["min_y"],
                                "max_x": extent["max_x"],
                                "max_y": extent["max_y"]
                            }
                        })

        return result