I have a CSV file with the Lat, Long and Rainfall Information. I would like to interpolate those point and create tiff file. Can any one can suggest me the easiest way to do that.
I am trying to using gdal_grid. I am very new on using gdal in python.
This is actually several questions. Assuming you have some scattered data for lats and longs you'll to build all the location were you want to make estimation (all lats and longs for the pixels of you Tiff image).
Once you have that you can use any of the solutions around to do IWD over your data (using a recent example in another question):
class Estimation():
# IWD. Check: https://stackoverflow.com/questions/36031338/interpolate-z-values-in-a-3d-surface-starting-from-an-irregular-set-of-points/36037288#36037288
def __init__(self,lon,lat,values):
self.x = lat
self.y = lon
self.v = values
def estimate(self,x,y,using='ISD'):
"""
Estimate point at coordinate x,y based on the input data for this
class.
"""
if using == 'ISD':
return self._isd(x,y)
def _isd(self,x,y):
#d = np.sqrt((x-self.x)**2+(y-self.y)**2)
d = x.copy()
for i in range(d.shape[0]):
d[i] = haversine(self.x[i],self.y[i],x,y)
if d.min() > 0:
v = np.sum(self.v*(1/d**2)/np.sum(1/d**2))
return v
else:
return self.v[d.argmin()]
The code above is actually adapted to calculate distance with the Haversine formula (which gives great-circle distances between two points on a sphere from their longitudes and latitudes). Notice again you can find all sorts of solutions for the haversine distance like this one:
def haversine(lon1, lat1, lon2, lat2):
"""
Check: https://stackoverflow.com/questions/15736995/how-can-i-quickly-estimate-the-distance-between-two-latitude-longitude-points
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
km = 6367 * c
return km
Finally once you have your array ready you should just build the Tiff using GDAL. For this check the following question for which I quote a part of it's solution:
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create('output.tif',xsize, ysize, 1, gdal.GDT_Float32, )
# this assumes the projection is Geographic lat/lon WGS 84
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
ds.SetProjection(srs.ExportToWkt())
gt = [ulx, xres, 0, uly, 0, yres ]
ds.SetGeoTransform(gt)
outband=ds.GetRasterBand(1)
outband.SetStatistics(np.min(mag_grid), np.max(mag_grid), np.average(mag_grid), np.std(mag_grid))
outband.WriteArray(mag_grid)
