Python: gridding point X,Y, and Z in order to extract statistical attribute

Sorry for the easy question but I am new in Python and i need same help.

My data are in point format: X,Y,Z. Where X and Y are coordinates, and z the value.

My problem is: create a raster (in TIF or ASCII) with 0.5 m by 0.5 m (or 1 by 1 m) where the value of each pixel is the avarage of Z. In case where i have not points in the pixel-i the value need to be NAN.

I am really glad for help with some code that i can study and implement,

Thanks in Advance for the help, I really need.

I tried to study and write a 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)

from this points i have problem to understand how index the x,y,z points in order to know where they drop. My first idea is give an index the the mesh grid and mark the points. After that I can make average for the points inside the pixel. Pixel empty (where there are not present points) are NAN.

But I don't know it's the right way to processing my data.

After that I wrote the following code to save in TIFF format 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

Really thanks for all help

Solution

A way to go:

define your grid spacing (a floating point number), which is the distance between two pixel/voxel midpoints in the same dimension
figure out the size of the grid you need, i.e. the number of grid points in x and y dimension, N_x, and N_y.
create two numpy arrays of that size with all values being zero using e.g. np.zeros([N_x, N_y])
iterate through your set of (x, y, v) points and
- project each (x, y) pair into its corresponding pixel, identified via two (integer) indices: x_i, y_i = tuple([int(c//spacing) for c in (x, y)])
- add a 1 to the one array at place (x_i, y_i) (holding the "count")
- add the value v to the other array at place (x_i, y_i) (holding the sum of values)
after having populated your two arrays, divide the sum-of-values-array by the count-array. 0/0.0 will automatically be assigned to NaN, while c/0.0 will be assigned to be Inf.