I have a set of points with X, and Y coordinates and associated Z values
POINT (-102.9071045 31.43453407 1871.930298)
POINT (-102.8339386 31.45111084 1856.14209)
POINT (-102.4687424 31.53144836 1869.777588)
POINT (-102.8522263 31.5140934 1852.520508)
POINT (-102.779068 31.53050041 1857.646606)
POINT (-102.7059631 31.54673767 1861.825195)
...... I am trying to create an interpolated surface of the Z values and save it as GeoTIFF. Here is the code I am using: '''
# Interpolate Point data to a raster grid
import os
import gdal
from osgeo import ogr
ogr.UseExceptions()
os.chdir('F:\PythonTutorial')
csv_file = 'Example.csv'
vrt_fn = 'Test12.vrt'
lyr_name = 'Example'
out_tiff = 'test14.tif'
#create a .vrt file
with open(vrt_fn, 'w') as fn_vrt:
fn_vrt.write('<OGRVRTDataSource>\n')
fn_vrt.write('\t<OGRVRTLayer name="%s">\n' % lyr_name)
fn_vrt.write('\t\t<SrcDataSource>%s.csv</SrcDataSource>\n' % lyr_name)
fn_vrt.write('\t\t<SrcLayer>%s</SrcLayer>\n' % lyr_name)
fn_vrt.write('\t\t<GeometryType>wkbPoint</GeometryType>\n')
fn_vrt.write('\t\t\t<LayerSRS>WGS84</LayerSRS>\n')
fn_vrt.write('\t\t<GeometryField encoding="PointFromColumns" srs="WGS84" x="Longitude" y="Latitude" z="CH4Concentration"/>\n')
fn_vrt.write('\t</OGRVRTLayer>\n')
fn_vrt.write('</OGRVRTDataSource>\n')
# print data on the screeen
inDataSource = ogr.Open(vrt_fn)
lyr = inDataSource.GetLayer(lyr_name)
for feat in lyr:
geom = feat.GetGeometryRef()
print(geom.ExportToWkt())
output = gdal.Grid(out_tiff, vrt_fn, algorithm = 'invdist:power=2.0:smoothing=1.0')
'''
The print data on the screen portion is added to check if the VRT file is being read - and it prints the point perfectly. However, the output image is null.
Can anyone help?