1

I'm using the world shapefile that can be found at https://trac.openstreetmap.org/browser/applications/editors/merkaartor-branches/merkaartor-0.14-fixes/share/world_shp/world_adm0.shp?rev=17052.

import pysal as ps
import ogr, osr
import shapely.geometry
from mpl_toolkits.basemap import Basemap

polygons = ps.open('world_adm0.shp')

# from http://gis.stackexchange.com/questions/129745/writing-a-shapefile-using-ogr-from-shapely-geometry-no-feature-added-error
# create the shapefile
driver = ogr.GetDriverByName("Esri Shapefile")
ds = driver.CreateDataSource("outputlocation.shp")
layr1 = ds.CreateLayer('',None, ogr.wkbPolygon)

for polygon in polygons:
    sh_poly = shapely.geometry.asShape(polygon)

    # create ogr geometry
    polylyr = ogr.CreateGeometryFromWkt(sh_poly.wkt)

    # create the field
    layr1.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))

    # Create the feature and set values
    defn = layr1.GetLayerDefn()
    feat = ogr.Feature(defn)
    feat.SetField('id', 1)
    feat.SetGeometry(polylyr)
    layr1.CreateFeature(feat)

# close the shapefile
ds.Destroy()

map = Basemap(width = 2000000, height = 1500000, resolution = 'l', projection = 'aea', lat_1 = -23, lat_2 = -32, lon_0 = 25, lat_0 = -29)
map.readshapefile('outputlocation', '', drawbounds=True)

results in:

ValueError                                Traceback (most recent call last)
<ipython-input-3-92ac4f1982b1> in <module>()
     33 
     34 map = Basemap(width = 2000000, height = 1500000, resolution = 'l', projection = 'aea', lat_1 = -23, lat_2 = -32, lon_0 = 25, lat_0 = -29)
---> 35 map.readshapefile('outputlocation', '', drawbounds=True)

/apps/python/2.7.6/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in readshapefile(self, shapefile, name, drawbounds, zorder, linewidth, color, antialiased, ax, default_encoding)
   2144         info = (shf.numRecords,shptype,bbox[0:2]+[0.,0.],bbox[2:]+[0.,0.])
   2145         npoly = 0
-> 2146         for shprec in shf.shapeRecords():
   2147             shp = shprec.shape; rec = shprec.record
   2148             npoly = npoly + 1

/apps/python/2.7.6/lib/python2.7/site-packages/mpl_toolkits/basemap/shapefile.pyc in shapeRecords(self)
    541         shapeRecords = []
    542         return [_ShapeRecord(shape=rec[0], record=rec[1]) \
--> 543                                 for rec in zip(self.shapes(), self.records())]
    544 
    545 class Writer:

/apps/python/2.7.6/lib/python2.7/site-packages/mpl_toolkits/basemap/shapefile.pyc in records(self)
    513         f.seek(self.__dbfHeaderLength())
    514         for i in range(self.numRecords):
--> 515             r = self.__record()
    516             if r:
    517                 records.append(r)

/apps/python/2.7.6/lib/python2.7/site-packages/mpl_toolkits/basemap/shapefile.pyc in __record(self)
    478                     value = float(value)
    479                 else:
--> 480                     value = int(value)
    481             elif typ == b('D'):
    482                 try:

ValueError: invalid literal for int() with base 10: '**********'

What exactly has gone wrong with writing the polygons?

1

2 Answers 2

2

You need the full shapefile, so better look into https://trac.openstreetmap.org/browser/subversion/applications/editors/merkaartor-branches/merkaartor-0.14-fixes/share/world_shp?rev=17052 and get all three files. Projection information is missing, so you have to assign WGS84 EPSG:4326 to it. Alternatively, take a full copy from Natural Earth.

The Albers equal area projection might fail on parts of Antarctica and along the dateline.

You have chosen the parameters for Africa (EPSG:102022 from ESRI). If that is your study area, clip the worldwide shapefile to that continent before reprojecting it to aea.

For worldwide coverage, you better choose another projection.

0

Is the shapefile EPSG:4326 ?

If not, you can convert via ogr2ogr -f "ESRI Shapefile" -t_srs "EPSG:4326" output.shp input.shp

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.