[GIS] Using PyShp to create polygon shapefiles

pyshp

I have a script for creating points shapefiles, but I can not fix it for make a new polyline shp file.

I have obtained the data from a excel file and extract data in three differents lists
Lon,lat,label

import shapefile as shp
from openpyxl import load_workbook
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
r = 'myfile.xlsm'


wb = load_workbook(filename= r, data_only=True, use_iterators = True)  
ws = wb['Export_Area']
cell_range = ws ['B5:D9']
tuple(ws.iter_rows('B5:D9'))
lista = []
for row in ws.iter_rows('B5:D9'):
    for cell in row:
        lista.append(cell.value)
        Longitud_area = (lista[1:len(lista):3])
        latitud_area = (lista[0:len(lista):3])
        puntos_area = (lista[2:len(lista):3])

creating the shp file

out_file = 'ExportArea.shp'


x,y,id_no=Longitud_area,latitud_area,puntos_area
logic = [True,False,True]
w = shp.Writer(shp.POINT)
w.autoBalance = 1 #ensures gemoetry and attributes match
w.field('X','F',10,5)
w.field('Y','F',10,5) #float - needed for coordinates
w.field('etiqueta','D') #date


for j,k in enumerate(x):
    w.point(k,y[j]) #write the geometry
    w.record(k,y[j],puntos_area[j]) #write the attributes

creating the prj file.

prj = open(out_file +'.prj', 'w')
proyeccion      ="GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'
prj.write(proyeccion)

It works fine for a point shp file but I can not see how to do the loop for making a poly shp file

Best Answer

If you want the coordinates (longitud, latitud) in a single list, use the function zip()

longitud= [-2.7692601473009155, -2.7362662697828974, -2.3913741501505776, -2.424387612708171, -2.7692601473009155]
latitud= [43.82769682706945, 43.66046724635432, 43.69596983959721, 43.86319942031234, 43.82769682706945]
coord = zip(longitud, latitud)
print coord
[(-2.7692601473009155, 43.82769682706945), (-2.7362662697828974, 43.66046724635432), (-2.3913741501505776, 43.69596983959721), (-2.424387612708171, 43.86319942031234), (-2.7692601473009155, 43.82769682706945)]

For reading an Excel file it is easiest now to use the Pandas module.

import pandas as pd
xl = pd.ExcelFile("mispuntos.xlsx")
print xl.sheet_names
[u'mispuntos']
# create a pandas dataframe
df = xl.parse('mispuntos')
print df.columns
Index([u'longitud', u'latitud', u'label'], dtype='object')
print df.head()
   longitud   latitud   label
0 -2.769260  43.827697     A
1 -2.736266  43.660467     B
2 -2.391374  43.695970     C
3 -2.424388  43.863199     D
4 -2.769260  43.827697   NaN

Point shapefile

import shapefile
w = shp.Writer(shp.POINT)
w.field('X','F',10,5)
w.field('Y','F',10,5) #float - needed for coordinates
w.field('label')
for index, row in df.iterrows():
   w.point(row['longitud'],row['latitud'])
   w.record(row['longitud'],row['latitud'],str(row['label']))
w.save('resulting')

Now, a polyline

# extract the 3 first rows of the dataframe
df2=df[0:3] # slicing
line = [[row['longitud'],row['latitud']] for index, row in df2.iterrows()]
print line # = list of points)
[[-2.7692601473, 43.8276968271], [-2.73626626978, 43.6604672464], [-2.39137415015, 43.6959698396]]
# write the polyline shapefile
w = shapefile.Writer(shapefile.POLYLINE)
w.field('label') 
w.line(parts=([line]))
w.record('a')
w.save('line')

Result

enter image description here

If you wish to go further, use GeoPandas (geospatial Pandas, the geometric operations are performed by shapely, file access by Fiona, not PyShp, and plotting by descartes and matplotlib)

Related Question