Geodjango：如何加载 .shp 文件并使用正确的 CRS 转换为 geojson？

Question

我有多个 shapefile (.shp) 及其辅助文件，我想将它们显示在 Leaflet 地图上。shapefile 使用不同的坐标参考系统 (CRS)，我很难掌握在地图上显示事物的最直接、最可靠的方式。在geodjango教程中，DataSource用于加载shapefile，然后对其进行操作。然而，在他们的示例中，他们只检索单个特征的几何形状，而不是整个 shapefile 的几何形状。我已经使用了PyShp，并且可以使用以下内容显示地图：

    sf = shapefile.Reader(filename)
    shapes = sf.shapes()
    geojson = shapes.__geo_interface__
    geojson = json.dumps(geojson)

但是，当 CRS 不是 WGS84 时，这会失败，事情不起作用，而且我不知道如何转换它。

阅读更多内容后，这篇文章抱怨了 CRS 支持和 pyshp，并建议使用 ogr2ogr。

因此，在尝试了解这些选项之后，我看到使用 Datasource、pyshp 和 ogr2ogr 作为可能的选项，但我不知道哪个选项真正最有意义。

我想要的只是使用 Django 将 .shp 文件转换为使用 WGS84 的 geojson 字符串，以便我可以将其包含在使用 Leaflet 的 HTML 页面上。

有经验的人可以推荐一条特定的路线吗？

Answer 1

没有一种直接的方法可以使用 Django 读取任何 shapefile DataSource，然后将其转换为EPSG:4326 (aka WGS84)，因此我们需要一步一步地创建并解决我们遇到的问题。

让我们开始这个过程：

1. .shp创建您需要读入的所有文件路径的列表。它应该如下所示：

   SHP_FILE_PATHS = [
       'full/path/to/shapefile_0.shp',
       'full/path/to/shapefile_1.shp',
       ...
       'full/path/to/shapefile_n.shp'
   ]
   

2. DataSource将 shapefile 读入对象。该信息存储在对象Layers（表示多层形状文件）中，该对象将其srs视为SpatialReference. 这很重要，因为我们稍后将转换几何图形以便WGS84可以在地图上显示。

3. 从每个 shapefile 的每一层，我们将利用该get_geoms()方法提取感知对象的列表。OGRGeometry srs

4. 每个这样的几何体都有一个json方法：

   以 JSON 格式返回此几何图形的字符串表示形式：

   >>> OGRGeometry('POINT(1 2)').json 
   '{ "type": "Point", "coordinates": [ 1.000000, 2.000000 ] }'
   

   这非常有用，因为这是创建FeatureCollection可在地图上显示的类型 geojson 的解决方案的一半。

5. geojsonFeatureCollection具有非常特定的格式，因此我们将创建基础并按程序填充它：

   feature_collection = {
       'type': 'FeatureCollection',
       'crs': {
           'type': 'name',
           'properties': {'name': 'EPSG:4326'}
       },
       'features': []
   }
   

6. 最后，我们需要features使用以下格式提取的几何图形填充列表：

   {
       'type': 'Feature',
       'geometry': {
           'type': Geometry_String,
           'coordinates': coord_list
        },
        'properties': {
            'name': feature_name_string
        }
   }
   

让我们将以上所有内容放在一起：

for shp_i, shp_path in enumerate(SHP_FILE_PATHS):
    ds = DataSource(shp_path)
    for n in range(ds.layer_count):
        layer = ds[n]
        # Transform the coordinates to epsg:4326
        features = map(lambda geom: geom.transform(4326, clone=True), layer.get_geoms())
        for feature_i, feature in enumerate(features):
            feature_collection['features'].append(
                {
                    'type': 'Feature',
                    'geometry': json.loads(feature.json),
                    'properties': {
                        'name': f'shapefile_{shp_i}_feature_{feature_i}'
                    }
                }
            )

现在feature_collection字典将包含转换为的提取的特征集合epsg:4326，您可以从它创建一个 json （例如json.dump(feature_collection)）

注意：虽然这可行，但似乎有点适得其反，您可以考虑将 shapefile 永久读入模型，而不是动态加载它们。