Python: point de maillage X, Y et Z pour extraire l'attribut statistique

Question

Désolé pour la question facile mais je suis nouveau en Python et j'ai besoin de la même aide. Mes données sont au format de point: X, Y, Z

Où X et Y sont des coordonnées et z la valeur. Mon problème est: créer un raster (en TIF ou ASCII) avec 0.5 m sur 0.5 m (ou 1 sur 1 m) où la valeur de chaque pixel est la moyenne de Z. Dans le cas où je n'ai pas de points dans le pixel-i la valeur doit être NAN.

Je suis vraiment heureux de l'aide avec un code que je peux étudier et mettre en œuvre,

Merci à l'avance pour l'aide, j'ai vraiment besoin.

J'ai essayé d'étudier et d'écrire un code:

from osgeo import gdal, osr, ogr 
import math, numpy 
import numpy as np 
import matplotlib.pyplot as plt 
import matplotlib.mlab as ml 
import matplotlib.delaunay 

tiff = 'test.tif' 

gridSize = 0.5 
# my area boundary 
xmax, xmin = 640000.06, 636999.83 
ymax, ymin = 6070000.3, 6066999.86 

# number of row and columns 
nx = int(math.ceil(abs(xmax - xmin)/gridSize)) 
ny = int(math.ceil(abs(ymax - ymin)/gridSize)) 

# Plot the points 
plt.scatter(x,y,c=z) 
plt.axis([xmin, xmax, ymin, ymax]) 
plt.colorbar() 
plt.show() 

# Generate a regular grid. 
xi = np.linspace(xmin, xmax, nx) 
yi = np.linspace(ymin, ymax, ny) 
xi, yi = np.meshgrid(xi, yi) 

de ce point que j'ai problème pour comprendre comment indexer le x, y, z des points afin de savoir où ils laissent tomber. Ma première idée est de donner un indice à la grille de maillage et de marquer les points. Après cela, je peux faire la moyenne pour les points à l'intérieur du pixel. Les pixels vides (où il n'y a pas de points présents) sont NAN.

Mais je ne sais pas si c'est la bonne façon de traiter mes données.

Après que j'ai écrit le code suivant pour enregistrer au format TIFF via GDAL

target_ds = gdal.GetDriverByName('GTiff').Create(tiff, nx,ny, 1, gdal.GDT_Byte) #gdal.GDT_Int32 
target_ds.SetGeoTransform((xmin, gridSize, 0,ymax, 0, -gridSize,)) 

if EPSG == None: 
    proj = osr.SpatialReference() 
    proj.ImportFromEPSG(EPSG) 
    # Make the target raster have the same projection as the source 
    target_ds.SetProjection(proj.ExportToWkt()) 
else: 
    # Source has no projection (needs GDAL >= 1.7.0 to work) 
    target_ds.SetProjection('LOCAL_CS["arbitrary"]') 

target_ds.GetRasterBand(1).WriteArray(numpy.zeros((ny,nx))) 
target_ds = None 

Vraiment merci pour toute l'aide

Answer 1

chemin à parcourir:

• définir votre grille spacing (un nombre à virgule flottante), qui est la distance entre deux points milieu pixel/voxel dans la même dimension
• taille de la grille dont vous avez besoin, c'est-à-dire le nombre de points de grille dans x et y dimension, N_x et N_y.
• Créez deux tableaux numpy de cette taille, toutes les valeurs étant nulles, par exemple. np.zeros([N_x, N_y])
• itérer sur votre ensemble de (x, y, v) des points et
  • projet chaque (x, y) paire dans son pixel correspondant, identifié par deux indices (entier): x_i, y_i = tuple([int(c//spacing) for c in (x, y)])
  • ajouter un 1 au tableau un à l'endroit (x_i, y_i) (tenant le « compte »)
  • ajouter la valeur v à l'autre tableau à l'endroit (x_i, y_i) (tenant la somme des valeurs)
• après avoir rempli vos deux tableaux , divise le tableau sum-of-values ​​par le count-array. 0/0.0 sera automatiquement affecté à NaN, tandis que c/0.0 sera affecté à Inf.