Convert gridded XYZ data (but without missing values) to GeoTIFF using gdal python interface. We needed to write this code because gdal_translate only works on xyz grid data when all the values of the grid, including missing values, are provided, and that's not the case here. Only the "good" values from the grid have been provided.
In [ ]:
# convert gridded XYZ data (but without missing values) to GeoTIFF
# data looks like this:
"""
379528.00 4577221.00 -5.38
379529.00 4577221.00 -5.38
379530.00 4577221.00 -5.40
379531.00 4577221.00 -5.41
379515.00 4577222.00 -5.34
379516.00 4577222.00 -5.36
"""
In [4]:
# UTM ZONE 19, NAD83 data
data =genfromtxt('/rps/bathy/muskegat/EW_Survey1.xyz')
In [5]:
x=data[:,0]
y=data[:,1]
z=data[:,2]
In [6]:
# find the extent and difference
xmin=min(x);xmax=max(x)
ymin=min(y);ymax=max(y)
mdx=abs(diff(x))
mdy=abs(diff(y))
In [7]:
# determine dx and dy from the median of all the non-zero difference values
dx=median(mdx[where(mdx>0.0)[0]])
dy=median(mdy[where(mdy>0.0)[0]])
In [8]:
#construct x,y,z of complete grid
xi=arange(xmin,xmax+dx,dx)
yi=arange(ymin,ymax+dy,dy)
zi=ones((len(yi),len(xi)))*NaN
shape(zi)
In [11]:
# calculate indices in full grid (zi) to stick the input z values
ix=numpy.round((x-xmin)/dx).astype(int)
iy=numpy.round((y-ymin)/dy).astype(int)
zi[iy,ix]=z
In [20]:
figure(figsize=(12,4))
# images start from the top left
zi=flipud(zi)
imshow(zi)
Out[20]:
In [22]:
# write as 32-bit GeoTIFF using GDAL
from osgeo import osr,gdal
ny,nx = zi.shape
driver = gdal.GetDriverByName("GTiff")
ds = driver.Create('output.tif', nx, ny, 1, gdal.GDT_Float32)
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
ds.SetGeoTransform( [ xmin, dx, 0, ymax, 0, dy ] )
# set the reference info
srs = osr.SpatialReference()
# UTM zone 19, North=1
srs.SetUTM(19,1)
srs.SetWellKnownGeogCS('NAD83')
ds.SetProjection( srs.ExportToWkt() )
# write the data to a single band
ds.GetRasterBand(1).WriteArray(zi)
# close
ds = None