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