Scripting en gvSIG: Análisis espacial (malla de puntos)

A falta de salir el acceso a imágenes raster desde el módulo de scripting, hemos hecho un pequeño script que simula esta función y permite realizar un análisis espacial o geoprocesamiento de la zona delimitada. Se basa en crear una malla de datos (elegiremos su densidad), y simular una capa de densidades.

En este caso, tenemos una simulación de una zona residencial muy sencilla, con zona de residencia, ocio, naturaleza, hospitales y carretera. Imaginaremos que queremos ver que zonas de residencia nos gustaría más para vivir, en función de las prioridades que tengamos respecto a las instalaciones, las podremos ajustar con unos porcentajes.

Estos rates vendrán en función de la distancia de cada punto a la instalación más cercana, quedándonos con el valor menor. Los rates podrían complicarse todo lo que queramos, desde hacerse en forma de puntuación, o si están en su área de influencia, o siguiendo una función, o tener en cuenta todas las geometrías u otras capas, etc..

Consiste en dos scripts, uno crea toda la malla y datos, y otro script sencillo recorta los que estén contenidos dentro de geometrías etiquetadas como residencia.

Si queréis probar el script que funciona directamente con la capa base1.shp podéis pedírmela por email que no me deja subirla a wordpress.

Si quieres adaptarlo a tu uso tendrás que ir buscando y cambiando sobretodo las lineas referidas al TIPO de campo de la geometría (ya sea natura, ocio, carretera…). En mi script aparece como: feature.campo1

crearMalla

from gvsig import *
from geom import *
from commonsdialog import *

def main():
    "Crea malla de puntos y calcula ratios de distancias a geometrias"
    #Capa base: base1.shp
    #campo1 con valores de:
    #"hospital","ocio","natura","carretera" y geometria

    layer = currentLayer()
    disPuntos = inputbox("Distancia entre puntos", "", QUESTION )
    disPuntos = int(disPuntos)
    if disPuntos <= 0:
        msgbox("El numero de puntos no es valido","AVISO",WARNING)
        return

    #Extensión de la capa
    envelope = layer.getFullEnvelope()
    xMin = int(envelope.getLowerCorner().getX())
    yMin = int(envelope.getLowerCorner().getY())
    xMax = int(envelope.getUpperCorner().getX())
    yMax = int(envelope.getUpperCorner().getY())
    values = dict()
    print "Creando malla.."

    #Datos para crear la capa
    ruta = "C:/gvsig/mdtarti/mallaRatios1.shp"
    CRS = layer.getProjectionCode()
    schema = createSchema()
    schema.append("GEOMETRY", "GEOMETRY")
    schema.append("ID", "INTEGER", size = 7)
    schema.append("X","INTEGER",size = 20)
    schema.append("Y","INTEGER",size = 20)
    schema.append("Valor", "DOUBLE", size = 7)
    schema.append("RatioSalud", "DOUBLE",size=10,default = 0)
    schema.append("RatioOcio", "DOUBLE",size=10,default = 0)
    schema.append("RatioNatura", "DOUBLE", size=10, default = 0)
    schema.append("RatioRoad", "DOUBLE", size = 10, default = 0)

    #Creamos la capa de salida
    output = createShape(schema,ruta,CRS=CRS,geometryType=POINT)

    #Calculamos la extension de la malla y puntos necesarios
    features= layer.features()
    rx = int(xMin)
    ry = int(yMin)
    index = 1
    numberPointsX = ((xMax - xMin) // disPuntos ) + 1
    numberPointsY = ((yMax - yMin) // disPuntos) + 1
    puntosTotales = numberPointsX*numberPointsY
    print "Puntos totales: ", puntosTotales
    #Restriccion de calculo maximo de puntos
    if puntosTotales > 10000: return
    if puntosTotales < 100: return

    #Calculo de la malla y de los ratios de cada punto
    for indiceY in range(numberPointsY):
        for indiceX in range(numberPointsX):
            #Posicion del punto
            values["X"] = rx
            values["Y"] = ry
            point = createPoint(int(values["X"]),int(values["Y"]))
            rx = rx + disPuntos

            #Calculo de los valores y ratios
            values["ID"] = index
            values["RatioSalud"] = ratio("hospital",point)
            values["RatioOcio"] = ratio("ocio",point)
            values["RatioNatura"] = ratio("natura",point)
            values["RatioRoad"] = ratio("carretera",point)
            values["GEOMETRY"] = point
            values["Valor"] = ratioValor(values["RatioSalud"],
                              values["RatioOcio"],
                              values["RatioNatura"],
                              values["RatioRoad"])
            index = index + 1

            #Añadimos la entidad
            output.append(values)
        ry = ry + disPuntos
        rx = int(xMin)
    print "Generado."
    output.commit()
    currentView().addLayer(output)

#Calculo del ratio segun la distancia
def ratio(campo,point):
      features = currentLayer().features()
      #Distancia maxima de influencia
      ratio = 10000
      for feature in features:
            if feature.campo1 == campo:
              #Si el punto esta dentro de la geometria
              if feature.geometry().contains(point):
                  ratio = 0
                  continue
              #Segun la distancia a la geometria
              if point.distance(feature.geometry()) < ratio:
                  ratio = point.distance(feature.geometry())
      return ratio

#Suma de ratios, personalizado segun lo que prefiramos
def ratioValor(salud,ocio,natura,carretera):
  porSalud = 0.20
  porOcio = 0.20
  porNatura = 1
  porCarretera = 0.30
  total = (salud * porSalud + ocio * porOcio +
    natura * porNatura + carretera * porCarretera)
  return total