I want to interpolate temperature in correlation with height. I have temperature data from stations with coordinates and height in this format:
Lat Lon Temp Height 46.2956 16.3978 10.71196 120 45.9776 16.0159 11.78539 150 45.8521 15.9598 11.18838 133 etc.
After reading the data there are 4 arrays (
height) and variable
z which represents correlation between temperature and height, e.g.,
z = -0.005 which means that for each meter of height temperature decreases for 0.005 °C.
Z is calculated like this:
z = np.polyfit(height, temp, 1,2)
After that I can calculate temperature on 0 meters for each station:
for i in range(len(temp)): temp0m.append((temp[i])-(height[i]*z))
I also have data with heights in meters with resolution of 0.003333333333 degrees and I want to use these heights and interpolated temperature to create map.
Heights are in
"dem_final2.txt" (ncolums = 2700, nrows = 1560)
1653 1571 1493 1429 1354... 1730 1699 1620 1528 1399...
So this is what I've done so far, everything is working but it's very slow:
if __name__ == '__main__': nx, ny = 1080,624 chunkSize = 10 t =  grid2 =  xEast = 12 ySouth = 42 xWest = 21 yNorth = 47.2 res = 0.003333333333 x,y,temp0m,z = downloadData() xi = np.linspace(xEast+0.1, xWest-0.1, nx) yi = np.linspace(ySouth+0.1, yNorth-0.1, ny) xi, yi = np.meshgrid(xi, yi) xi, yi = xi.flatten(), yi.flatten() l = len(xi)/chunkSize for i in range(chunkSize): interp = rbfInterp(x,y,temp0m,xi[(l*i):(l*(i+1))],yi[(l*i):(l*(i+1))]) grid2 = grid2 + list(interpData(xi[(l*i):(l*(i+1))],yi[(l*i):(l*(i+1))],z,interp,xEast,ySouth,xWest,yNorth,res))[(l*i):(l*(i+1))] grid2 = np.array(grid2, np.float) grid2 = grid2.reshape((ny, nx)) plot (x,y,temp0m,grid2) #plotting to png with matplotlib def rbfInterp(x, y, temp0m,xi,yi): interp = Rbf(x, y, temp0m, function='linear') return interp(xi,yi) def interpData(xi,yi,z,interp,xEast,ySouth,xWest,yNorth,res): append = t.append row= -1 for i in xrange(len(xi)): tmp = int(round((yi[i]-ySouth)/res,0)) tmp1 = int(round((xi[i]-xEast)/res,0)+1) if (row!=tmp): tmp2 = linecache.getline('dem_final2.txt', tmp).split(' ') if (float(tmp2[tmp1])>-10): append ((float(tmp2[tmp1]))*z+float(interp[i])) else: append (None) row = tmp return t
I divided into chunks because a lot of memory was used without them and I think this can be used to speed up things maybe with more processes or something like that but I don't have knowledge to do that without any help. Of course, any other help is also welcome.
Btw this script runs for about 25-30 sec