5

I've got some scattered data in the form of (latitude, longitude, someParameterValue). I'm using inverse distance weighting interpolation method to interpolate them in a rectangular grid of pixels. Presently I'm generating the query points for that grid, in python, as given below. Please note that I've converted the (latitude, longitude) coordinates to cartesian (x, y) coordinates :

xr = int(math.ceil(xmax-xmin)); yr = int(math.ceil(ymax-ymin)) xres=yres=1 xr = math.ceil(xr/xres)+1 yr = math.ceil(yr/yres)+1 npts = int(xr * yr) queryPts = np.zeros(shape=(npts,2)) x1=xmin; y1=ymin idx=0 while(x1 <= xmax): while(y1 <= ymax): #the 2D array of query points is populated here queryPts[idx] = [x1, y1] idx += 1 y1 += yres y1=ymin x1 += xres

where xmin, ymin, xmax, ymax are the minimum and maximum values of x and y coordinates respectively. Here I feel that populating the query points at intervals of 1 in each of x and y axes is not the right way to go. After calculating the grid of interpolated values, I'm using gdal to turn it into a raster image with the interpolated values scaled to 0-255 for the pixels. A sample image, that I assume to be inappropriate, is shown here. interpolated image for the query points calculated using the above step

Would like to get suggestions on the following:

  1. What would be the proper way to generate query points for an interpolation grid i.e. how to set the resolution of the points on x and y axes whose interpolated values that we are going to calculate?
  2. I can see that the above question is related to the pixel resolution that we want in the final interpolated image. Hence the above question could be asked as: how to generate 2D arrays of (x,y) points (query points?) that correspond to pixels in the rendered image?
2

i may try to attempt the answers in reverse:

the resolution might be somewhat specific to what you want - and as you mention, dependent on the cell size. If your coordinates were in meters - you might want to end with 10 meter cell size, so you might use

cellsize = 10
ncol = int(math.ceil(xmax-xmin)) / cellsize 
nrow = int(math.ceil(ymax-ymin)) / cellsize
xres = (xmax - xmin) / float(ncol) #which should get back to original cell size
yres = (ymax - ymin) / float(nrow)

to generate a 2D array of points for interpolation, i sometimes use numpy mgrid as

import numpy as np
gridx, gridy = np.mgrid[xmin:xmax:ncol*1j, ymin:ymax,nrow*1j]

note the step argument is a complex number in the example above - from the docs 'However, if the step length is a complex number (e.g. 5j), then the integer part of its magnitude is interpreted as specifying the number of points to create between the start and stop values, where the stop value is inclusive.' of course, you could also use linspace, or mesh grid, or probably any number of functions to generate the x and y points.

i found this post to be very helpful. and used it to evaluate different methods of interpolation - where 'pts' is a list of x,y pairs and z is a list of z values - ptx and pty are lists of x and y values respectively.

interptype = 'gauss': #or rbf or griddata

 ##### using griddata #####
 if interptype == 'griddata':
     from scipy.interpolate import griddata
     grid = griddata(pts,z,(gridx,gridy), method='linear',fill_value=-3e30)

 ##### using radial basis function ####
 if interptype == 'rbf':
     import scipy.interpolate as interpolate
     f = interpolate.Rbf(ptx, pty, z, function='linear')
     grid = f(gridy, gridx)

 ##### using gaussian ####
 if self.interptype == 'gauss':
     from sklearn.gaussian_process import GaussianProcess
     ptx = np.array(ptx)
     pty = np.array(pty)
     z = np.array(z)
     print math.sqrt(np.var(z))
     gp = GaussianProcess(regr='quadratic',corr='cubic',theta0=np.min(z),thetaL=min(z),thetaU=max(z),nugget=0.05)
     gp.fit(X=np.column_stack([pty,ptx]),y=z)
     rr_cc_as_cols = np.column_stack([gridy.flatten(), gridx.flatten()])
     grid = gp.predict(rr_cc_as_cols).reshape((ncol,nrow))

i think the resolution doesn't show up again until writing the raster (something like this?) -

outTiff = 'somefilename.tif'
geotransform=(xmin,xres,0,ymax,0, -yres)   
# That's (top left x, w-e pixel resolution, rotation (0 if North is up), 
#         top left y, rotation (0 if North is up), n-s pixel resolution)
drv = gdal.GetDriverByName('GTiff')
ds = drv.Create(outTiff,ncol, nrow, 1 ,gdal.GDT_Float32)  # Open the file
band = ds.GetRasterBand(1)
band.SetNoDataValue(-3e30)
ds.SetGeoTransform(geotransform)  # Specify its coordinates
if wkt != '':
     ds.SetProjection(wkt)   # Exports the coordinate system 
band.WriteArray(np.flipud(grid.T)) # if you need to flip or transpose your grid
  • Thank you @fluidmotion. After seeing mgrid, I tried different interpolation methods. It works well! I'm looking to plot this interpolated map on top of Google maps, any suggestion with respect to KML layers (or any standard procedure for that) would be great. If needed, I'll post that as a separate question. – freax Jun 14 '15 at 11:40
  • i haven't exported raster to kml (at least not in a long time) - a quick search turns up some potentially helpful sites - gdal2tiles may be the way to go. and an older post asking something similar.. – fluidmotion Jun 14 '15 at 13:28
  • Thanks! I see that the whole interpolation and raster image conversion process is slow (around 1-2 seconds in a system with RAM 8G for approx 5000 scattered points). Given that I'd like to run this interpolation process on a browser, what would be best solution? Javascript look more attractive, at least if some one were to animate the generated images but for that I'd have to find apt counter-parts of the calculations that are done in python. – freax Jun 29 '15 at 11:06
  • not sure about the best browser method - might be best as a new question to get more suggestions? – fluidmotion Jun 29 '15 at 12:10

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.