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?

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.