Introducción

Si quieres leer o escribir un archivo shp en Python, probablemente vayas a usar pyshp (o puede que Fiona).

He intentado encontrar un complemento de Python para la libreria pyshp en pypi con el que pudiera realizae queries en la base de datos contenida en el archivo .dbf.

Dado la sección censal española obtenemos el shapefile (INE Cartografía digitalizada) y lo visualizamos con Mapshaper:

El objetivo es seleccionar algunas partes del mapa del shapefile y salvarlo como otra colección de archivos .shp, .dbf, shx. Por ejemplo:

select from Spanish census section:
- province: Madrid
- autonomy: Castilla y León, Galicia

Shapefile

Shapefile es un formato de datos vectoriales geoespaciales. Realmente, es una colección de tres archivos:

`.shp: geometria binaria (poligonos…).
.dbf: datos de las figuras o registros. En formato dBase.
.shx: índice de figuras (no es obligatorio) para indexar más rápido.

También es posible que se incluyan otros archivos como .proj, .shp.xml, .sbn …

Query snippet

El archivo .dbf tiene los siguientes campos para cada figura:

['OBJECTID', 'CUSEC', 'CUMUN', 'CSEC', 'CDIS', 'CMUN', 'CPRO', 'CCA', 'CUDIS', 'OBS', 'CNUT0', 'CNUT1', 'CNUT2', 'CNUT3', 'CLAU2', 'NPRO', 'NCA', 'NMUN', 'Shape_Leng', 'Shape_area', 'Shape_len']

Por ejemplo, NPROV es el nombre de la provincia, NMUN es el nombre del municipio yCUMUN es el código del municipio.

Queremos selección en la misma query múltiples campos con diversos valores. Por ejemplo, NPRO = Madrid, Sevilla y NMUN = Barcelona.

import shapefile

class ShapeFileUtils(shapefile.Reader):
    ''' Inherit from shapefile.Reader '''
    
    def __init__(self,  *args, **kwargs):
        shapefile.Reader.__init__(self,  *args, **kwargs)
        # list of name fields
        self.name_fields = [field[0] for field in
                            list(self.fields[1:])]
    
    def query(self, **args):
        '''
        Makes a query on a shapefile.Reader object
        :param args: list of fields and values to query
        :return: the query shapefile.Writer object
        '''
        w = shapefile.Writer(shapeType=self.shapeType)
        w.fields = self.fields[1:]
        # shapefile 2.x version:
        # w.encoding = self.encoding
        for key, value in args.items():
            index = self.name_fields.index(key)
            for shaperec in self.iterShapeRecords():
                if shaperec.record[index] in value:
                    w.record(*shaperec.record)
                    # shapefile 2.x version should use:
                    # w.shape(shaperec.shape)
                    # shapefile 1.x version should use:
                    w._shapes.append(shaperec.shape)
        return w

Primero, cargamos el shapefile del [Instituto Nacional de Estadísitca]](http://www.ine.es/censos2011_datos/cen11_datos_resultados_seccen.htm) español en el objeto ShapeFileItilsque extiende de la clase shapefile.Reader:

sf = ShapeFileUtils("SECC_CPV_E_20111101_01_R_INE",
                     encoding="latin1")

Después, realizamos una query seleccionamos algunas provincias del Norte de España, la Comunidad Autónoma de Castilla y León y algunos de los municipios más poblados de Asturias, pero sin seleccionar la provincia:

wq = sf.query(
         NPRO=['Coruña, A', 'Ourense', 'Lugo', 'Pontevedra', 'Cantabria'],
         NMUN=['Oviedo', 'Gijón', 'Avilés'],
         NCA=['Castilla y León']
    )
<shapefile.Writer object at 0x111278908>
wq.save('shapefiles/test/query_test')

El codigo anterior devolverá un objeto shapefile.Writer que podremos visualizar en Mapshaper:

También es posible realizar dos quieries para guardar dos shapefiles disintos y representarlos juntos. En este caso el municipio de Madrid y su provincia:

wq2 = sf.query(NPRO=['Madrid'])
wq3 = sf.query(NMUN=['Madrid'])

Vista de la sección censal del Retiro, dentro del municipio y la provincia de Madrid:

Próximos pasos

Resulta algo incómodo realizar las queris de esta manera. Estaría bien poder realizar consultas de manera simiar a cómo se hace con el ORM de Django. Por ejemplo:

sf.objects.filter( NPRO=’Barcelona’ ).exclude( NMUN=’Prat de Llobregat, El’ )…