Your approach is good, but you could make things clearer using dictionaries instead of list and the csv module allows it.
Moreover, your script use two loops while it is possible to simplify using only one (the second loop is redundant, one line of the csv file = one record of the shapefile).
1) With dictionaries:
Reading the csv file:
with open('your.csv', 'rb') as f:
reader = csv.DictReader(f)
for row in reader:
print row
{'xr': '22.5', 'yb': '31.353634', 'name': 'some Name', 'xl': '-25.3125', 'yt': '47.517193'}
{'xr': '-0.703125', 'yb': '74.40216', 'name': 'another Name', 'xl': '-103.359375', 'yt': '80.87282'}
You can now use row['xr'] or row['name'] instead of row[n] (more explicit)
So your script becomes:
import csv
polyName, polyPart = [],[]
with open('your.csv', 'rb') as f:
reader = csv.DictReader(f)
for row in reader:
bl = [float(row['xl']),float(row['yb'])]
tl = [float(row['xl']),float(row['yt'])]
br = [float(row['xr']),float(row['yb'])]
tr = [float(row['xr']),float(row['yt'])]
parr = [tl, tr, br, bl, tl]
polyName.append(row['name'])
polyPart.append(parr)
Writing the polygon shapefile
If you look at PyShpDocs you can see that:
A polygon is defined by:
w = shapefile.Writer(shapefile.POLYGON)
w.line(parts=[[[1,5],[5,5],[5,1],[3,3],[1,1]]])
and the script, as you wrote is:
import shapefile
# create the Polygon shapefile
w = shapefile.Writer(shapefile.POLYGON)
# the field
w.field('name','C',maxStringLength)
# write the polygons in the shapefile
for part,name in zip(polyPart, polyName):
w.poly(parts=[part])
w.record(name)
#save the shapefile
w.save('your.shp')
2) Final solution using only one loop
But the second loop is not necessary here: you can do it all with one loop (reading the csv file and writing the shapefile without using the polyName and polyPart lists).
w = shapefile.Writer(shapefile.POLYGON)
w.field('name','C',50)
with open('your.csv', 'rb') as f:
reader = csv.DictReader(f)
for row in reader:
bl = [float(row['xl']),float(row['yb'])]
tl = [float(row['xl']),float(row['yt'])]
br = [float(row['xr']),float(row['yb'])]
tr = [float(row['xr']),float(row['yt'])]
parr = [tl, tr, br, bl, tl]
w.poly(parts=[parr])
w.record(row['name'])
w.save("your.shp')
Result in QGIS:
The correct solution is to use m.proj4stringand not directly albers_crs because it is not a "valid" Proj4 string (even if it apparently works: major axis or radius = 0 not given)
print m.proj4string
'+lon_0=10.0 +y_0=1250000.0 +R=6370997.0 +proj=aea +x_0=2000000.0 +units=m +lat_2=55.0 +lat_1=35.0 +lat_0=45.0 '
# so
pyproj_convert = Proj(m.proj4string)
pyproj_x, pyproj_y = pyproj_convert(lon, lat)
print 'pyproj:', pyproj_x, pyproj_y
pyproj: 1932284.3542 1653858.27802
You can control the result with the Python module GDAL
from osgeo import osr
from_srs = osr.SpatialReference()
from_srs.ImportFromEPSG(4326)
to_srs = osr.SpatialReference()
to_srs.ImportFromProj4(m.proj4string)
transf = osr.CoordinateTransformation(from_srs,to_srs)
print transf.TransformPoint(lon,lat)[:2]
(1932284.3541970258, 1653858.278018097)
Best Answer
Use directly the GDAL module rather than using it through (Geo)Django:
or the pyproj module (with
preserve_units=Truebecause pyproj assumes that your coordinates are in meters and this is not the case for EPSG:2227):Combine with Fiona: