1

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?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.