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

Question

我有多个 shapefile (.shp) 及其要显示在 Leaflet map 上的辅助文件。 shapefile 使用不同的坐标引用系统 (CRS)，我很难掌握在 map 上显示事物的最直接和可靠的方式。在 geodjango tutorial , DataSource 用于加载一个shapefile，然后对其进行操作。然而，在他们的例子中，他们只检索单个特征的几何，而不是整个 shapefile。我用过 PyShp我可以使用以下内容显示 map :

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

但是，当 CRS 不是 WGS84 时，这将失败，并且我不知道如何转换它。
多读一点， this post提示 CRS 支持和 pyshp，并建议使用 ogr2ogr。
所以在尝试了解选项后，我看到使用 Datasource、pyshp 和 ogr2ogr 作为可能的选项，但我不知道哪个选项最有意义。
我想要的只是将使用 Django 的 .shp 文件转换为使用 WGS84 的 geojson 字符串，以便我可以将其包含在使用 Leaflet 的 HTML 页面中。
任何有更多经验的人都可以推荐一条特定的路线吗？

Answer 1

没有使用 Django 的 DataSource 读取任何 shapefile 的直接方法。然后将其翻译为 EPSG:4326 (aka WGS84) ，因此我们需要一步一步地创建并解决我们遇到的问题。
让我们开始这个过程:

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

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

DataSource将 shapefile 读入一个对象。该信息存储在对象的 Layers 中(代表一个多层 shapefile)知道他们的 srs作为 SpatialReference . 这很重要，因为我们稍后会将几何转换为 WGS84以便在 map 上显示。

从每个 shapefile 的每一层，我们将利用 get_geoms() 提取 OGRGeometry 列表的方法srs有意识的对象。

每个这样的几何体都有一个 json 方法:

Returns a string representation of this geometry in JSON format:

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

这非常有用，因为它是创建 FeatureCollection 的解决方案的一半输入将显示在 map 上的 geojson。

一个 FeatureCollection geojson 有一个非常特殊的格式，因此我们将创建基础并按程序填充它:

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

最后，我们需要填充 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 dict 将包含转换为 epsg:4326 的提取特征集合您可以创建一个 json 表单(例如 json.dump(feature_collection) )
注意:虽然这会起作用，但似乎有点适得其反，您可以考虑将 shapefile 永久读入模型，而不是即时加载它们。