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?

nmtoken
  • 13,355
  • 5
  • 38
  • 87
johnsmith
  • 141
  • 4
  • Did you download just the .shp file, you also need the .dbf file and the .shx file from https://trac.openstreetmap.org/browser/subversion/applications/editors/merkaartor-branches/merkaartor-0.14-fixes/share/world_shp?rev=17052 – nmtoken Nov 01 '15 at 10:29

2 Answers2

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.

AndreJ
  • 76,698
  • 5
  • 86
  • 162
0

Is the shapefile EPSG:4326 ?

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

sam-6174
  • 111
  • 4