GeoDataFrame多边形EPSG:32748转EPSG:4326报ProjError求助

问题背景 #

有一个名为poi_shp_utm_buffered的GeoDataFrame，存储多边形数据，CRS为EPSG:32748。调用poi_shp_utm_buffered.to_crs(4326)转换至EPSG:4326时触发上述错误。

GeoDataFrame数据结构 #

|   |                                   PlaceID | Official_N | Primary__1 | Full_Stree | Overall_Sc |             X |            Y |                                          geometry |
|--:|------------------------------------------:|-----------:|-----------:|-----------:|-----------:|--------------:|-------------:|--------------------------------------------------:|
| 0 | 360qqu5b-409576faa57505ae78aa1cd551661af6 |          a |          1 |          1 |          2 | 709469.120296 | 9.299854e+06 | POLYGON ((709479.120 9299853.811, 709479.072 9... |
| 1 | 360qqu5b-5ec43f6ad613e60c15e3c4779f9c003d |          b |          1 |          1 |          2 | 709369.905462 | 9.299615e+06 | POLYGON ((709379.905 9299615.157, 709379.857 9... |
| 2 | 360qqu5b-7c11918bf9754eb6841c838f3337b783 |          c |          2 |          2 |          1 | 707546.465569 | 9.300030e+06 | POLYGON ((707556.466 9300030.011, 707556.417 9... |

CRS信息 #

>>> poi_shp_utm_buffered.crs

<Projected CRS: EPSG:32748>
Name: WGS 84 / UTM zone 48S
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: World - S hemisphere - 102°E to 108°E - by country
- bounds: (102.0, -80.0, 108.0, 0.0)
Coordinate Operation:
- name: UTM zone 48S
- method: Transverse Mercator
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

报错回溯 #

执行代码：

poi_shp_utm_buffered.to_crs(4326)

报错信息：

ProjError                                 Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_15024\2103948270.py in 
      2 # __ = poi_shp_utm_buffered
      3 
----> 4 __ = poi_shp_utm_buffered.drop(columns = ['X','Y']).to_crs(4326)

c:\Users\arasyidi\Anaconda3\envs\python_ds_gis\lib\site-packages\geopandas\geodataframe.py in to_crs(self, crs, epsg, inplace)
   1362         else:
   1363             df = self.copy()
-> 1364         geom = df.geometry.to_crs(crs=crs, epsg=epsg)
   1365         df.geometry = geom
   1366         if not inplace:

c:\Users\arasyidi\Anaconda3\envs\python_ds_gis\lib\site-packages\geopandas\geoseries.py in to_crs(self, crs, epsg)
   1122         """
   1123         return GeoSeries(
-> 1124             self.values.to_crs(crs=crs, epsg=epsg), index=self.index, name=self.name
   1125         )
   1126 

c:\Users\arasyidi\Anaconda3\envs\python_ds_gis\lib\site-packages\geopandas\array.py in to_crs(self, crs, epsg)
    777         transformer = Transformer.from_crs(self.crs, crs, always_xy=True)
    778 
-> 779         new_data = vectorized.transform(self.data, transformer.transform)
...
    432             iny,

pyproj/_transformer.pyx in pyproj._transformer._Transformer._transform()

ProjError: x, y, z, and time must be same size

环境版本信息 #

>>> gpd.show_versions()

SYSTEM INFO
-----------
python     : 3.9.13 (main, Aug 25 2022, 23:51:50) [MSC v.1916 64 bit (AMD64)]
executable : c:\Users\arasyidi\Anaconda3\envs\python_ds_gis\python.exe
machine    : Windows-10-10.0.22000-SP0

GEOS, GDAL, PROJ INFO
--------------------- 
GEOS       : None
GEOS lib   : None
GDAL       : 3.6.2
GDAL data dir: c:\Users\arasyidi\Anaconda3\envs\python_ds_gis\lib\site-packages\pyogrio\gdal_data\
PROJ       : 6.2.1
PROJ data dir: C:\Users\arasyidi\Anaconda3\envs\python_ds_gis\Library\share\proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.12.2
numpy      : 1.21.5
pandas     : 1.4.3
pyproj     : 2.6.1.post1
shapely    : 1.8.4
fiona      : None
geoalchemy2: None
geopy      : None
matplotlib : 3.5.2
mapclassify: 2.5.0
pygeos     : 0.9
pyogrio    : 0.5.1
psycopg2   : None
pyarrow    : 11.0.0
rtree      : 0.9.7

补充信息 #

使用conda list查看环境时显示已安装fiona，但gpd.show_versions()无法检测到它。

解决方法 #

1. 检查并清理几何数据
   该错误多因Geometry列存在空值、无效几何或维度不匹配要素，执行以下代码排查：

   # 检查空几何
   print(poi_shp_utm_buffered[poi_shp_utm_buffered.geometry.isna()])
   # 检查无效几何
   print(poi_shp_utm_buffered[~poi_shp_utm_buffered.geometry.is_valid])
   

   清理异常数据：

   poi_clean = poi_shp_utm_buffered.dropna(subset=['geometry']).copy()
   poi_clean = poi_clean[poi_clean.geometry.is_valid]
   

2. 更新依赖库版本
   当前pyproj版本（2.6.1.post1）较旧，与GDAL 3.6.2存在兼容性风险，执行更新：

   conda update pyproj geopandas
   

3. 指定pyogrio处理转换
   由于fiona未被检测到，强制使用pyogrio：

   from pyogrio import set_default_driver
   set_default_driver('GPKG')
   poi_wgs84 = poi_shp_utm_buffered.to_crs(epsg=4326)
   

4. 手动实现坐标转换
   若上述方法无效，用pyproj手动转换：

   from pyproj import Transformer
   import shapely
   
   transformer = Transformer.from_crs("EPSG:32748", "EPSG:4326", always_xy=True)
   
   def transform_geom(geom):
       if geom is None or not geom.is_valid:
           return None
       coords = list(transformer.itransform(geom.exterior.coords))
       return shapely.Polygon(coords)
   
   poi_shp_utm_buffered['geometry'] = poi_shp_utm_buffered.geometry.apply(transform_geom)
   

内容的提问来源于stack exchange，提问作者Amri Rasyidi