根据经纬度获取当前所在城市————自己开发城市逆地址解析2(城市数据导入方法)

上篇文章讲述了根据经纬度坐标，查询所在城市的接口实现方法。本篇讲述，如何将下载的城市geojson格式的轮廓，导入到postgresql数据库中。

示例开源地址:zzgeoapi

导入遇到的问题

1、geojson中城市中心“centrer”类型为OFTRealList，在geodjango的LayerMapping内无法对应PointField，“bbox”也希望有类型无法对应的问题。

解决方法:自定义GeojsonLayerMapping，在导入时使用轮廓自己计算轮廓中心和bbox(参见_set_extern_fields方法)

2、有些城市轮廓有些问题，在判断内点时会报错。

解决方法:引入shapely库，在导入数据时判断轮廓是否可用，如果有问题，使用shapely.make_valid进行修复，修复后的GeometryCollection存入mpoly2中。

废话不说，上代码

import os
import sys
from pathlib import Path
from django.contrib.gis.utils import LayerMapping, LayerMapError
from django.contrib.gis.gdal import DataSource
from django.core.exceptions import FieldDoesNotExist, ObjectDoesNotExist
from django.contrib.gis.db.models import GeometryField
from django.contrib.gis.gdal import (
    CoordTransform, DataSource, GDALException, OGRGeometry, OGRGeomType,
    SpatialReference,
)
from django.db import connections, models, router, transaction
from django.contrib.gis.db.models.functions import MakeValid
from django.contrib.gis.geos.error import GEOSException
# from django.contrib.gis.geos import Polygon, MultiPolygon, GeometryCollection

import shapely


from apps.zzgeo.models import AdArea

# Build paths inside the project like this: BASE_DIR / 'subdir'.
BASE_DIR = Path(__file__).resolve().parent.parent


ERR_GEOS = []


class GeojsonLayerMapping(LayerMapping):
    def check_layer(self):
        """
        Check the Layer metadata and ensure that it's compatible with the
        mapping information and model. Unlike previous revisions, there is no
        need to increment through each feature in the Layer.
        """
        # The geometry field of the model is set here.
        # TODO: Support more than one geometry field / model.  However, this
        # depends on the GDAL Driver in use.
        self.geom_field = False
        self.fields = {}

        # Getting lists of the field names and the field types available in
        # the OGR Layer.
        ogr_fields = self.layer.fields
        ogr_field_types = self.layer.field_types

        # Function for determining if the OGR mapping field is in the Layer.
        def check_ogr_fld(ogr_map_fld):
            try:
                idx = ogr_fields.index(ogr_map_fld)
            except ValueError:
                raise LayerMapError('Given mapping OGR field "%s" not found in OGR Layer.' % ogr_map_fld)
            return idx

        # No need to increment through each feature in the model, simply check
        # the Layer metadata against what was given in the mapping dictionary.
        for field_name, ogr_name in self.mapping.items():
            # Ensuring that a corresponding field exists in the model
            # for the given field name in the mapping.
            try:
                model_field = self.model._meta.get_field(field_name)
            except FieldDoesNotExist:
                raise LayerMapError('Given mapping field "%s" not in given Model fields.' % field_name)

            # Getting the string name for the Django field class (e.g., 'PointField').
            fld_name = model_field.__class__.__name__

            if isinstance(model_field, GeometryField):
                if self.geom_field:
                    raise LayerMapError('LayerMapping does not support more than one GeometryField per model.')

                # Getting the coordinate dimension of the geometry field.
                coord_dim = model_field.dim

                try:
                    print('ogr_name: ' + ogr_name)
                    if coord_dim == 3:
                        gtype = OGRGeomType(ogr_name + '25D')
                    else:
                        gtype = OGRGeomType(ogr_name)
                except GDALException as ee:
                    raise LayerMapError('Invalid mapping for GeometryField "%s".' % field_name)

                # Making sure that the OGR Layer's Geometry is compatible.
                # ltype = self.layer.geom_type
                # if not (ltype.name.startswith(gtype.name) or self.make_multi(ltype, model_field)):
                #     raise LayerMapError('Invalid mapping geometry; model has %s%s, '
                #                         'layer geometry type is %s.' %
                #                         (fld_name, '(dim=3)' if coord_dim == 3 else '', ltype))

                # Setting the `geom_field` attribute w/the name of the model field
                # that is a Geometry.  Also setting the coordinate dimension
                # attribute.
                self.geom_field = field_name
                self.coord_dim = coord_dim
                fields_val = model_field
            elif isinstance(model_field, models.ForeignKey):
                if isinstance(ogr_name, dict):
                    # Is every given related model mapping field in the Layer?
                    rel_model = model_field.remote_field.model
                    for rel_name, ogr_field in ogr_name.items():
                        idx = check_ogr_fld(ogr_field)
                        try:
                            rel_model._meta.get_field(rel_name)
                        except FieldDoesNotExist:
                            raise LayerMapError('ForeignKey mapping field "%s" not in %s fields.' %
                                                (rel_name, rel_model.__class__.__name__))
                    fields_val = rel_model
                else:
                    raise TypeError('ForeignKey mapping must be of dictionary type.')
            else:
                # Is the model field type supported by LayerMapping?
                if model_field.__class__ not in self.FIELD_TYPES:
                    raise LayerMapError('Django field type "%s" has no OGR mapping (yet).' % fld_name)

                # Is the OGR field in the Layer?
                idx = check_ogr_fld(ogr_name)
                ogr_field = ogr_field_types[idx]

                # Can the OGR field type be mapped to the Django field type?
                if not issubclass(ogr_field, self.FIELD_TYPES[model_field.__class__]):
                    raise LayerMapError('OGR field "%s" (of type %s) cannot be mapped to Django %s.' %
                                        (ogr_field, ogr_field.__name__, fld_name))
                fields_val = model_field

            self.fields[field_name] = fields_val
            self.layer_code = int(self.layer.name)

    def _valid_mpoly(self, m):
        try:
            m.mpoly.contains(m.mpoly.centroid)
        except GEOSException:
            ERR_GEOS.append(m.mpoly)
            print('=================geo contains error==================')
            for poly in m.mpoly:
                for ring in poly:
                    print(list(ring))
            print('=================geo contains error==================')
            import shapely
            polys = []
            for poly in m.mpoly:
                _p = [shapely.Polygon(item) for item in poly]
                polys.extend(_p)
            valid_mpoly = shapely.make_valid(shapely.MultiPolygon(polys))
            from shapely.geometry.collection import GeometryCollection as _GeometryCollection
            from shapely.geometry import MultiPolygon as _MultiPolygon
            if isinstance(valid_mpoly, shapely.geometry.collection.GeometryCollection):
                setattr(m, 'mpoly2', OGRGeometry(valid_mpoly.wkt).wkt)
                m.mpoly2.contains(m.mpoly.centroid)
            elif isinstance(valid_mpoly, _MultiPolygon):
                setattr(m, 'mpoly', OGRGeometry(valid_mpoly.wkt).wkt)
                m.mpoly.contains(m.mpoly.centroid)
            else:
                print('valid_mpoly class: %s' % valid_mpoly.__class__)
                raise

    def _set_extern_fields(self, m):
        self._valid_mpoly(m)
        m.parent_code = self.layer_code
        m.center = m.mpoly.centroid
        m.bbox = m.mpoly.envelope

        print(m.parent_code)
        print(m.center)
        print(m.bbox)

    def save(self, verbose=False, fid_range=False, step=False,
             progress=False, silent=False, stream=sys.stdout, strict=False):
        default_range = self.check_fid_range(fid_range)
        # Setting the progress interval, if requested.
        if progress:
            if progress is True or not isinstance(progress, int):
                progress_interval = 1000
            else:
                progress_interval = progress

        def _save(feat_range=default_range, num_feat=0, num_saved=0):
            if feat_range:
                layer_iter = self.layer[feat_range]
            else:
                layer_iter = self.layer

            for feat in layer_iter:
                code = feat.get('code')
                if code == 0 or code is None:
                    continue
                num_feat += 1
                # Getting the keyword arguments
                try:
                    kwargs = self.feature_kwargs(feat)
                except LayerMapError as msg:
                    # Something borked the validation
                    if strict:
                        raise
                    elif not silent:
                        stream.write('Ignoring Feature ID %s because: %s\n' % (feat.fid, msg))
                else:
                    # Constructing the model using the keyword args
                    is_update = False
                    if self.unique:
                        # If we want unique models on a particular field, handle the
                        # geometry appropriately.
                        try:
                            # Getting the keyword arguments and retrieving
                            # the unique model.
                            u_kwargs = self.unique_kwargs(kwargs)
                            m = self.model.objects.using(self.using).get(**u_kwargs)
                            is_update = True

                            # Getting the geometry (in OGR form), creating
                            # one from the kwargs WKT, adding in additional
                            # geometries, and update the attribute with the
                            # just-updated geometry WKT.
                            geom_value = getattr(m, self.geom_field)
                            if geom_value is None:
                                geom = OGRGeometry(kwargs[self.geom_field])
                            else:
                                geom = geom_value.ogr
                                new = OGRGeometry(kwargs[self.geom_field])
                                for g in new:
                                    geom.add(g)
                            setattr(m, self.geom_field, geom.wkt)
                        except ObjectDoesNotExist:
                            # No unique model exists yet, create.
                            m = self.model(**kwargs)
                    else:
                        print(type(kwargs.get('mply')))
                        m = self.model(**kwargs)

                    try:
                        # Attempting to save.
                        self._set_extern_fields(m)
                        m.save(using=self.using)
                        num_saved += 1
                        if verbose:
                            stream.write('%s: %s\n' % ('Updated' if is_update else 'Saved', m))
                    except Exception as msg:
                        if strict:
                            # Bailing out if the `strict` keyword is set.
                            if not silent:
                                stream.write(
                                    'Failed to save the feature (id: %s) into the '
                                    'model with the keyword arguments:\n' % feat.fid
                                )
                                stream.write('%s\n' % kwargs)
                            raise
                        elif not silent:
                            stream.write('Failed to save %s:\n %s\nContinuing\n' % (kwargs, msg))

                # Printing progress information, if requested.
                if progress and num_feat % progress_interval == 0:
                    stream.write('Processed %d features, saved %d ...\n' % (num_feat, num_saved))

            # Only used for status output purposes -- incremental saving uses the
            # values returned here.
            return num_saved, num_feat

        if self.transaction_decorator is not None:
            _save = self.transaction_decorator(_save)

        nfeat = self.layer.num_feat
        if step and isinstance(step, int) and step < nfeat:
            # Incremental saving is requested at the given interval (step)
            if default_range:
                raise LayerMapError('The `step` keyword may not be used in conjunction with the `fid_range` keyword.')
            beg, num_feat, num_saved = (0, 0, 0)
            indices = range(step, nfeat, step)
            n_i = len(indices)

            for i, end in enumerate(indices):
                # Constructing the slice to use for this step; the last slice is
                # special (e.g, [100:] instead of [90:100]).
                if i + 1 == n_i:
                    step_slice = slice(beg, None)
                else:
                    step_slice = slice(beg, end)

                try:
                    num_feat, num_saved = _save(step_slice, num_feat, num_saved)
                    beg = end
                except Exception:  # Deliberately catch everything
                    stream.write('%s\nFailed to save slice: %s\n' % ('=-' * 20, step_slice))
                    raise
        else:
            # Otherwise, just calling the previously defined _save() function.
            _save()



area_mapping = {
    'code': 'code',
    'name': 'name',
    'fullname': 'fullname',
    # 'center': 'center',
    'children_num': 'childrenNum',
    'level': 'level',
    # 'bbox': 'bbox',
    "mpoly": "MULTIPOLYGON",
}


def import_ds(ds, verbose):
    print(ds)
    lm = GeojsonLayerMapping(AdArea, ds, area_mapping, transform=False)
    lm.save(strict=True, verbose=verbose)



def run(verbose=True):
    china_geojson = os.path.join(BASE_DIR, 'data', 'geojson', '100000.json')
    ds = DataSource(china_geojson)
    import_ds(ds, verbose)
    for item in ds[0]:
        code = item.get('code')
        if code == 0 or code is None:
            continue
        item_geojson = os.path.join(BASE_DIR, 'data', 'geojson', '%d.json' % code)
        if not os.path.exists(item_geojson):
            continue
        item_ds = DataSource(item_geojson)
        import_ds(item_ds, verbose)

导入方法

1、运行django环境

python manage.py shell

2、导入脚本

from tools import load_china_city_data as load

load.run()

 

本文抛砖引玉，欢迎讨论