数据分析或处理中关于坐标系的一些事

 2011年，GOCE卫星的一幅新大地水准面地图。

想象中的地球

常用地理坐标系：

投影坐标系 (PCS) 是使用线性测量而非角度单位的平面系统。投影坐标系由地理坐标系和地图投影共同组成。地图投影包含数学计算，用于将 GCS 的角度大地坐标转换为平面 PCS 系统的笛卡尔坐标。

投影坐标系 = 地理坐标系 + 投影方法

常用投影：

（1）

这种圆锥投影使用两条标准纬线，相比仅使用一条标准纬线的投影，可在某种程度上减少畸变。此投影最适用于中纬度东西方向分布的大陆板块。

阿尔伯斯投影由 Heinrich C. Albers 于 1805 年提出。椭圆体方程由 Oscar S.Adams 于 1927 年开发。

当你使用pyproj库中的CRS对象时，可以直接通过EPSG代码来初始化这个对象，例如：

from pyproj import CRS# 创建一个基于EPSG代码的CRS对象crs = CRS.from_epsg(3857)  # 这是Web Mercator投影的EPSG代码

这样，pyproj就利用EPSG代码加载了对应的坐标参考系统定义，并且可以使用这个定义来进行坐标转换操作。例如，你可以使用pyproj.Transformer来创建一个转换器，将坐标从一个EPSG系统转换到另一个EPSG系统：

from pyproj import Transformer# 创建一个转换器，用于从EPSG:4326（WGS 84）到EPSG:3857（Web Mercator）的转换transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857")lon, lat = -117.0, 34.0x, y = transformer.transform(lon, lat)

@classmethoddef calculate_epsg(cls, longitude, coords_code, zone_width):    cls._check_params(longitude, coords_code, zone_width)    if coords_code == 4490:        base_epsg = 4509 if zone_width == 6 else 4548        base_central_meridian = 117    elif coords_code == 4326:        base_epsg = None  # This case does not apply to the EPSG calculation logic used.        base_central_meridian = None    elif coords_code == 4610:        base_epsg = 2343 if zone_width == 6 else 2380        base_central_meridian = 105    elif coords_code == 4214:        base_epsg = 21458 if zone_width == 6 else 2432        base_central_meridian = 105    elif coords_code == 4555:        base_epsg = 4584 if zone_width == 6 else 4792        base_central_meridian = 105    else:        base_epsg = None        base_central_meridian = None    if base_epsg is not None and base_central_meridian is not None:        central_meridian = cls._calculate_zone(longitude, zone_width)        epsg_code = cls._calculate_epsg(base_epsg, base_central_meridian, central_meridian, zone_width)    elif coords_code == 4326:        pp = int(longitude / 6) + 31        epsg_code = f"326{pp}"    else:        # epsg_code = 102100        epsg_code = 3395    return epsg_code

@classmethoddef create_custom_projection(cls, center_longitude, coords_code=4490, zone_width=1.5):    """Create a custom projection based on the provided parameters."""    center_longitude = cls._calculate_zone(center_longitude, zone_width)    # Retrieve ellipsoid parameters    if coords_code in code_params:        params = code_params[coords_code]        a = params["a"]        rf = params["rf"]    else:        print(f"Ellipsoid '{coords_code}' not found in ellipsoid_params.")        return None    # Construct the projection definition string    '''    tmerc: 横向墨卡托投影    merc: 墨卡托投影（正轴）    lcc: 兰勃脱圆锥投影（Lambert Conformal Conic）    sterea: Stereographic projection    utm: UTM投影（实际上是基于横向墨卡托，但有一些特殊的参数设置）    '''    center_longitude = int(float(center_longitude)) if float(center_longitude).is_integer() else float(center_longitude)    ellipsoid = code_name[coords_code]    custom_proj_definition = (        f"+proj=tmerc +lat_0=0 +lon_0={center_longitude} +k=1 +x_0=500000 +y_0=0 "        f"+a={a} +rf={rf} +units=m +no_defs +proj_info:({ellipsoid} / {zone_width}-degree Gauss-Kruger CM {center_longitude}E) +coords_code:[{coords_code}]"    )    try:        custom_crs = CRS(custom_proj_definition)    except Exception as e:        print(f"Error creating CRS: {e}")        return None    return custom_crs

def cal_distance_between_line(line1, line2, spa_type):    if spa_type == "Projected":        return line1.distance(line2)    else:        geom1, geom2 = nearest_points(line1, line2)         # 近似计算测地线        return geodesic((geom1.y, geom1.x), (geom2.y, geom2.x)).meters

def cal_distance(geom1, geom2, spatial_ref):    spa_type = spatial_ref.type    if spa_type == "Projected":        # 将 Point 对象转换为 PointGeometry 对象，然后计算距离        point_geom1 = arcpy.PointGeometry(geom1, spatial_ref)        point_geom2 = arcpy.PointGeometry(geom2, spatial_ref)        return point_geom1.distanceTo(point_geom2)    else:        return geodesic((geom1.Y, geom1.X), (geom2.Y, geom2.X)).meters

class Geographic2Project(GeoFC, FormatEpsg):    def __init__(self, _in_fc, zone_width):        super().__init__(_in_fc)        self.fc_name = os.path.basename(self.fc).split('.')[0]        self.gdf = self.read_fc(include_oid=False)        self.zone_width = zone_width        self.gdf_code = self.gdf.crs.to_epsg()        self.geo_type = self.gdf.geometry.geom_type.unique()[0].lower()        self.cal_group()    # 按照 epsg_code 分类进行投影到投影坐标系    @classmethod    def apply_projection(cls, group, transformer, crs):        # 应用投影函数到整个几何列        group['geometry'] = group['geometry'].apply(cls.project_geometry, args=(transformer, cls.geo_type))        # 处理投影后的几何体        group['geometry'] = group['geometry'].apply(cls.process_geometry)        # 移除空几何体        group = group.dropna(subset=['geometry'])        # 生成一个新的 gdf，包含投影坐标系中的几何信息        projected_gdf = gpd.GeoDataFrame(group, geometry='geometry', crs=crs)        return projected_gdf if not projected_gdf.empty else None    def cal_group(self):        # 确认 gdf 的坐标系类型        if self.gdf.crs.is_projected:            raise ValueError("GeoDataFrame is already in a projected coordinate system. No projection needed.")        elif self.gdf.crs.is_geographic:            gdf_code = self.gdf.crs.to_epsg()            geo_type = self.gdf.geometry.geom_type.unique()[0].lower()            # 计算每个点位的投影带 EPSG 代码            if geo_type == "point":                self.gdf['longitude'] = self.gdf.geometry.x            else:                self.gdf['longitude'] = self.gdf.geometry.apply(lambda geom: geom.centroid.x)            if self.zone_width in [3, 6]:                self.gdf['epsg_code'] = self.gdf['longitude'].apply(self.calculate_epsg, coords_code=gdf_code,                                                                    zone_width=self.zone_width)            else:                self.gdf['epsg_code'] = self.gdf['longitude'].apply(                    lambda lon: str(self._calculate_zone(lon, zone_width=self.zone_width)))        else:            raise ValueError("Unsupported CRS type. Only projected or geographic coordinate systems are supported.")    def geographic_to_project(self):        # 按照 epsg_code 分类进行投影到投影坐标系        projected_gdfs = {}        for epsg_code, group in self.gdf.groupby('epsg_code'):            # 解析epsg_code，用于区分按3°、6°，还是其他经度跨度进行的投影            if abs(float(epsg_code)) > 180:                epsg_code = int(epsg_code)                custom_proj = f"EPSG:{epsg_code}"                transformer = Transformer.from_proj(f"EPSG:{self.gdf_code}", custom_proj, always_xy=True)            else:                center_longitude = int(float(epsg_code)) if float(epsg_code).is_integer() else float(epsg_code)                custom_proj = self.create_custom_projection(center_longitude, self.gdf_code, zone_width=self.zone_width)                transformer = Transformer.from_proj(f"EPSG:{self.gdf_code}", custom_proj, always_xy=True)            # 按照 epsg_code 分类进行投影到投影坐标系            projected_gdf = self.apply_projection(group, transformer, custom_proj)            projected_gdfs[epsg_code] = projected_gdf        # 移除空值（跳过 EPSG:4326）        projected_gdfs = {k: v for k, v in projected_gdfs.items() if v is not None}        return projected_gdfs    def project_export(self, out_path=None, out_name=None, columns_to_remove=None, export=False):        projected_gdfs = self.geographic_to_project()        if export is True:            if out_path is None:                raise ValueError("输出目录不能为空！")            export_obj = ExportData()            export_obj.gdf_to_gdb(projected_gdfs, out_path, out_name, columns_to_remove, fc_name=self.fc_name)        return projected_gdfs

class Project2Geographic(FormatEpsg):    def __init__(self, gdf_or_dict):        super().__init__()        self.gdf_or_dict = gdf_or_dict    def project_to_geographic(self):        if isinstance(self.gdf_or_dict, dict):            # Create a defaultdict to store GeoDataFrames grouped by geographic CRS code            grouped_gdfs = defaultdict(list)            # Iterate over each key-value pair in the dictionary            for zone, gdf in self.gdf_or_dict.items():                if gdf.crs.is_projected:                    target_code = self.get_geographic_crs_code(gdf.crs)                    gdf = gdf.to_crs(target_code)                    grouped_gdfs[target_code].append(gdf)                elif gdf.crs.is_geographic:                    grouped_gdfs[gdf.crs.to_epsg()].append(gdf)                else:                    raise ValueError(                        f"Unsupported CRS type for zone {zone}. Only projected or geographic coordinate systems are supported.")            # Process each group of GeoDataFrames            result_gdfs = {}            for target_code, group_gdfs in grouped_gdfs.items():                # Combine all GeoDataFrames in the same group into one                combined_gdf = gpd.GeoDataFrame(pd.concat(group_gdfs, ignore_index=True), crs=f'EPSG:{target_code}')                result_gdfs[int(target_code)] = combined_gdf            return result_gdfs        elif isinstance(self.gdf_or_dict, gpd.GeoDataFrame):            gdf = self.gdf_or_dict            if gdf.crs.is_projected:                target_code = self.get_geographic_crs_code(gdf.crs)                geo_type = gdf.geometry.geom_type.unique()[0].lower()                transformer = Transformer.from_crs(gdf.crs.to_epsg(), f"EPSG:{target_code}", always_xy=True)                if geo_type == "point":                    def project_geometry(geometry):                        lon, lat = transformer.transform(geometry.x, geometry.y)                        return Point(lon, lat)                else:                    def project_geometry(geometry):                        project = partial(transformer.transform, direction='forward')                        return transform(project, geometry)                gdf['geometry'] = gdf['geometry'].apply(project_geometry)                gdf.crs = f'EPSG:{target_code}'                return gdf            elif gdf.crs.is_geographic:                return gdf            else:                raise ValueError("Unsupported CRS type. Only projected or geographic coordinate systems are supported.")        else:            raise ValueError("Input must be either a GeoDataFrame or a dictionary of GeoDataFrames.")    def geographic_export(self, out_path, out_name=None, columns_to_remove=None, export=False):        geographic_gdfs = self.project_to_geographic()        if export is True:            if out_path is None:                raise ValueError("输出目录不能为空！")            export_obj = ExportData()            export_obj.gdf_to_gdb(geographic_gdfs, out_path, out_name, columns_to_remove)        return geographic_gdfs

if __name__ == '__main__':    in_feature = r"D:\7.python 工具箱\工具箱部署\3.测试数据\水系线\云南省_水系线数据_4490.shp"    zone_width = 3    out_path = r"D:\7.python 工具箱\工具箱部署\3.测试数据\水系线\scratch.gdb"    out_name = None    columns_to_remove = ['longitude', 'epsg_code']    # 对输入要素按投影带进行分组，分组投影后导出    to_project_obj = Geographic2Project(in_feature, zone_width)    project_gdfs = to_project_obj.project_export(out_path, out_name, columns_to_remove, export=True)    # 又将分组投影后导出后的结果合并，变换为地理坐标系导出    to_geographic_obj = Project2Geographic(project_gdfs)    project_gdfs = to_geographic_obj.geographic_export(out_path, out_name, columns_to_remove=None, export=True)

扫码关注我们

从进阶开始，共同探索

GIS高级玩法……

扫描左侧二维码加微信

QQ：376512832

添加微信或QQ请备注：

GIS探案