Suite

Quel pourcentage du polygone x est dans le pixel raster y, efficacement


Je pense que c'est un problème digne d'une compétition, en raison de sa taille.

Je voudrais pondérer un ensemble de pixels raster par rapport aux attributs d'un ensemble de polygones. Pour ce faire, j'ai besoin de connaître la proportion de chaque polygone dans chaque pixel raster.

J'aimerais également conserver les identifiants du pixel et du polygone sur la sortie finale.

Parfois les polygones sont plus grands que les pixels, parfois plus petits. Les pixels avec une valeur de 0 peuvent être exclus.

Les polygones avec lesquels je travaille sont Statistical Area Level 1 (SA1) ASGS Ed 2011 et ils peuvent être trouvés au bas de cette page :

http://www.abs.gov.au/AUSSTATS/[email protected]/DetailsPage/1270.0.55.001July%202011?OpenDocument

Le raster avec lequel je travaille est ici :

http://www.abs.gov.au/AUSSTATS/[email protected]/DetailsPage/1270.0.55.0072011?OpenDocument

Je crois que ce sont GDA94, les polygones sont epsg:4283, le raster est epsg:3577

Comme vous le remarquerez, ils sont tous deux volumineux, c'est pourquoi l'efficacité est importante.

Toute solution open source (R, QGIS, etc.) serait la meilleure, bien que s'il existe une solution efficace dans ArcGIS ou MapInfo, il serait intéressant de la connaître.


d'après ma compréhension de votre question et en regardant les données, vous souhaitez attribuer une population à chaque polygone. En pseudocode, quelque chose comme

pour chaque polygone popn <- 0 pour chaque intersection (polygone, contour de pixel) fraction <- calculer la zone d'intersection (en pourcentage de la valeur du contour de pixel) popn <- popn + (population*fraction)

Pour ce faire, vous pouvez utiliser QGIS et PostgreSQL/PostGIS

Dans QGIS,

  • convertir le raster en polygones. Assurez-vous de cocher la case « DN », afin que les valeurs de pixels soient copiées. Notez que ce ne sera pas un mappage 1:1 entre les pixels et les polygones - les pixels contigus avec la même valeur formeront une zone. Plus à ce sujet plus tard…

  • sélectionnez tous les polygones avec une valeur > 0 et enregistrez la sélection en tant que nouveau fichier de formes. Cela accélérera les choses en supprimant la plupart des Aus, qui sont en grande partie vides ;)

En zoomant, vous devriez maintenant voir quelque chose comme ça

Notez que les formes irrégulières ont une valeur unique - par ex. 3. Cela signifie que chaque pixel de la forme a une valeur de 3, donc si la forme a une surface équivalente à 5 pixels, cette forme a en fait une population de 15.

  • importer ce shapefile dans postgres en utilisant shp2pgsql, en tant que table 'pixels'

  • faire de même pour les zones statistiques, comme le tableau "zones"

J'ai utilisé epsg:3577, tout ce qui compte, c'est que les deux sont identiques.

Vous pouvez maintenant utiliser postgis pour obtenir l'intersection de chaque polygone de zone avec chaque polygone de "zone de pixel" et résumer les valeurs de chaque polygone de

(a/b) * (pop*(b/p))
  • où a est l'aire de l'intersection avec l'aire du pixel
  • où b est la zone de la zone de pixel (qui peut-être un pixel, ou plusieurs)
  • où p est l'aire d'un seul pixel (constante, 1 km^2)
  • où pop est la valeur (population) de la zone de pixels

qui devrait simplifier jusqu'à

pop * (a/p)

Alors, quelque chose comme

créer un chevauchement de table en tant que ( sélectionnez a.sa1_main11, p.DN * (ST_Area(ST_Intersection(a.geom,p.geom))/1000000.0) en tant que weighted_pop, ST_Intersection(a.geom,p.geom) à partir de pixels en tant que p, zones en tant que where a.geom && p.geom limit 10000 );

En ramenant cela dans QGIS et en coloriant chaque intersection au hasard, vous pouvez voir comment les zones ont été morcelées.

Vous devriez pouvoir voir que les populations totalisent les mêmes valeurs que l'original.

Vous pouvez désormais effectuer une requête pour agréger les populations pondérées sur le champ sa1_main11. Vous pouvez joindre votre fichier de formes de zone d'origine au résultat de cette requête.

évidemment, enlevez la limite lorsque vous voulez tout faire. Je l'ai fait pour pouvoir obtenir un résultat en quelques minutes à vérifier… mais il s'est produit dans la majeure partie de l'Australie occidentale.


Juste au cas où quelqu'un serait venu ici pour un problème comme celui-ci, voici quelques notes :

• epsg : 3577 (comme l'a utilisé Steven Kay) est la référence spatiale correcte dans ce cas (Australie) car elle préserve l'aire. Les références spatiales qui ne préservent pas exactement la zone (comme epsg:4283) renverront une zone légèrement erronée. Tout devra être re-projeté sur la zone choisie en préservant la projection (merci aaryno de l'avoir signalé.)

• Les pixels contigus avec la même valeur ne doivent pas former une zone - car bien que la population non pondérée soit la même, la population pondérée sera différente.

• Je pense que les fonctions correctes dans PostGIS sont la géographie (sphérique) plutôt que la géométrie (cartésienne).

• L'avantage d'utiliser la géométrie est qu'elle est beaucoup moins complexe et donc plus rapide à exécuter - au point où c'est faisable. Au moins sous les contraintes de mon système, la géographie n'est pas faisable, même limitée à 10 000, sans parler des 500 000 impairs requis.

• La solution que j'ai utilisée consiste à calculer le chevauchement en fonction de la géométrie non projetée selon la réponse de Steven Kay. Ensuite, dans un deuxième temps, reprojetez les polygones résultants en « géographie » et calculez la zone.