Geodjango: how to load .shp file and convert to geojson with the right CRS?

user_1162541 0
Viewed 1785

I have multiple shapefiles (.shp) with their auxiliary files that I want to display on a Leaflet map. The shapefiles use different coordinate reference systems (CRS) and I struggle to grasp the most straightforward and reliable way to show things on the map. In the geodjango tutorial, DataSource is used to load a shapefile and then manipulate it. However, in their examples, they only retrieve the geometry of individual features, not of the entire shapefile. I have used PyShp and I am able to show a map using something like:

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

However, this fails when the CRS is not WGS84 things don't work, and I don't see how to convert it.

Reading a bit more, this post complains about CRS support and pyshp, and suggests using ogr2ogr.

So after trying to understand the options, I see using Datasource, pyshp, and ogr2ogr as possible options, but I don't know which option really makes the most sense.

All I want is convert a .shp file using Django into a geojson string that uses WGS84 so that I can include it on an HTML page that uses Leaflet.

Can anyone with more experience suggest a particular route?

geodjango coordinate-transformation shapefile django python
2 Answers

There is not a straight forward way to read any shapefile using Django's DataSource and then translate it to EPSG:4326 (aka WGS84), therefore we need to create one step by step and solve the issues that arise as we get to them.

Let's begin the process:

  1. Create a list of all the .shp file paths that you need to read in. That should look like this:

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

  2. DataSource reads the shapefile into an object. The information is stored in the object's Layers (representing a multilayered shapefile) that are aware of their srs as a SpatialReference. That is important because we will transform the geometry later to WGS84 in order to be displayable on the map.

  3. From each Layer of each shapefile, we will utilize the get_geoms() method to extract a list of OGRGeometry srs aware objects.

  4. Each such geometry has a json method that:

    Returns a string representation of this geometry in JSON format:

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

    That is very useful because is half the solution to create a FeatureCollection type geojson that will be displayable to the map.

  5. A FeatureCollection geojson has a very specific format and therefore we will create the basis and fill it procedurally:

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

  6. Finally, we need to populate the features list with the extracted geometries in the following format:

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

Let's put all the above together:

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}'
                    }
                }
            )

Now the feature_collection dict will contain the extracted feature collection transformed to epsg:4326 and you can create a json form it (ex. json.dump(feature_collection))

NOTE: Although this will work, it seems a bit counterproductive and you may consider reading the shapefiles into a model permanently instead of loading them on the fly.

user_74149390 
#models.py

from __future__ import unicode_literals
from django.db import models
from django.contrib.gis.db import model
# county


class County(models.Model):
    district = models.CharField(max_length=50)
    count = models.FloatField()
    county_nam = models.CharField(max_length=50)
    code = models.BigIntegerField()
    geom = models.MultiPolygonField(srid=4326)

    def _unicode_(self):
        return self.county_nam

    class Meta:
        verbose_name_plural="County"

in views

from django.shortcuts import render
from django.views.generic import TemplateView
from django.core.serializers import serialize
from django.http import HttpResponse
from .models import County

# Create your views here.

class HomePageView(TemplateView):
    template_name='index.html'

# county views

def county_datasets(request):
    county=serialize('geojson',County.objects.all())
    return HttpResponse(county,content_type='application/json

in urls.py

from django.conf.urls import url
from djgeojson.views import GeoJSONLayerView
from .views import HomePageView,county_datasets

# create your urls here


urlpatterns=[
    url(r'^$',HomePageView.as_view(),name='home'),
    url(r'^county_data/$',county_datasets,name='county'),

]

run the server

python manage.py runserver

check your geojson data by typing

http://localhost:8000/county_data/

user_123590530
Related