0

ArcGIS does not employ a pure linear interpolation algorithm, so I turned to python and specifically scipy to perform the interpolation.

I found out I can create an ASCII file containing the interpolated values and pass that to ArcMap in order to come up with a DEM. I successfully performed the interpolation, then created the ASCII per the instruction provided in the ESRI ASCII Raster format documentation and converted the ASCII to a raster using ArcMap's corresponding tool.

I first tried importing random values in order to check that this procedure actually works and it worked as expected. But when I pass the interpolated values, although the conversion tool does not raise any errors, the DEM makes no sense.

Below is the python code I used and an image of the DEM.

Can anyone solve any issues or recommend a workaround or another way to do it?

import sys
import numpy as np
from scipy.interpolate import griddata


np.set_printoptions(formatter={'float_kind':'{:f}'.format})
x,y,z = np.loadtxt(sys.argv[1],delimiter=',',usecols=(1,2,3),unpack=True)
pixel = np.float(sys.argv[3])
xmin = min(x)
ymin = min(y)
xmax = max(x)
ymax = max(y)
nrows = (ymax - ymin)/pixel
ncols = (xmax - xmin)/pixel
dir_txt = sys.argv[1].replace('.txt', '')

xi,yi = np.mgrid[xmin:xmax:complex(0, ncols), ymin:ymax:complex(0, nrows)]
zi = griddata((x, y), z, (xi, yi), fill_value=-9999, method=sys.argv[2])

# ASCII file header
header = "ncols %d \nnrows %d\nxllcorner %f\nyllcorner %f\ncellsize %f\nnodata_value -9999" % (int(ncols), int(nrows), float(min(x)), float(min(y)), pixel)
np.savetxt("out.txt", zi, fmt='%f', header=header, comments='', delimiter=' ', newline='\n')

DEM created from ASCII file

Improve this question
5
  • Can you show us the actual values input to create your output? If there are too many values for any reason, then can you simplify your data to a smaller sample too, please?
    – PolyGeo ♦
    Commented Feb 3, 2018 at 22:13
  • The exact points that I used to interpolate can be found here and the output ASCII file can be found here.
    – Nikos
    Commented Feb 4, 2018 at 7:56
  • I think they should be included in your question (not as a link) - the output could then be generated by anyone volunteering their time to test.
    – PolyGeo ♦
    Commented Feb 4, 2018 at 8:08
  • @PolyGeo thank you for the vast amount of help so far. You helped me figure out a workaround by using the griddata function provided by matplotlib instead of scipy and that solved my problems. Thanks again for your invaluable help!
    – Nikos
    Commented Feb 4, 2018 at 8:43
  • I think it would still be helpful to the community if you were to include that input data in your question, and to add your solution as an answer.
    – PolyGeo ♦
    Commented Feb 4, 2018 at 8:46

2 Answers 2

Reset to default
1

With your second solution (with matplotlib)

enter image description here

But you can get the same result with Scipy (with the file you supplied) 

import pandas as pd
import numpy as np
test = pd.read_csv("sample.txt")
test.head()
    id        x            y       z
0  5000  304388.973  4239207.555  14.503
1  5001  304383.335  4239210.993  14.274
2  5002  304387.942  4239207.943  14.070
3  5003  304386.662  4239208.707  14.649
4  5004  304384.396  4239210.264  14.658
xmin,xmax,ymin,ymax = [min(test.x),max(test.x),min(test.y),max(test.y)]
pixel = 10
nrows = int((ymax - ymin)/pixel)
ncols = int((xmax - xmin)/pixel)

With your original solution 

xi,yi = np.mgrid[xmin:xmax:complex(0, ncols), ymin:ymax:complex(0, nrows)]
zi = il.griddata((test.x, test.y), test.z, (xi, yi),fill_value=-9999, method='linear')

enter image description here

The problem is due to the calculation of the grid (np.mgrid(). With a simple np.meshgrid(

xi = np.linspace(xmin, xmax, ncols)
yi = np.linspace(ymin, ymax, nrows)
xi, yi = np.meshgrid(xi, yi) 
zi = il.griddata((test.x, test.y), test.z, (xi, yi),fill_value=-9999, method='linear')

enter image description here

You need to flip the numpy array (with numpy.flipup)

zi = np.flipud(zi)
header = "ncols %d \nnrows %d\nxllcorner %f\nyllcorner %f\ncellsize %f\nnodata_value -9999" % (ncols,nrows, xmin, ymin, pixel)
np.savetxt("out5.asc", zi, fmt='%f', header=header, comments='', delimiter=' ', newline='\n')

enter image description here

Improve this answer
1
  • Thank you very much! I will accept your answer because I prefer working with SciPy. I came up with some triangulation issues when I used big datasets with matplotlib.
    – Nikos
    Commented Feb 5, 2018 at 9:10
0

After hours of trial and error as well as google searching, I came up with a solution. I used the corresponding matplotlib griddata function instead of the one provided by SciPy. Below is a working code, a point sample in PENZ format and the output ASCII file.

import sys
import numpy as np
import matplotlib.mlab as ml


np.set_printoptions(formatter={'float_kind':'{:f}'.format})
x,y,z = np.loadtxt(sys.argv[1],delimiter=',',usecols=(1,2,3),unpack=True)
pixel = np.float(sys.argv[2])
xmin = min(x)
ymin = min(y)
xmax = max(x)
ymax = max(y)
nrows = (ymax - ymin)/pixel
ncols = (xmax - xmin)/pixel

dir_txt = sys.argv[1].replace('.txt', '_interp.asc')

xi = np.linspace(min(x), max(x), num=ncols)
yi = np.linspace(min(y), max(y), num=nrows)
zi = ml.griddata(x, y, z, xi, yi, interp='linear')

zi = np.rot90(np.fliplr(zi), 2)

zi[np.isnan(zi)]=-9999

# ASCII file header
header = "ncols %d \nnrows %d\nxllcorner %f\nyllcorner %f\ncellsize %f\nnodata_value -9999" % (int(ncols), int(nrows), float(min(x)), float(min(y)), pixel)
np.savetxt(dir_txt, zi, fmt='%f', header=header, comments='', delimiter=' ', newline='\n')

Sample Point File:

5000,304388.973,4239207.555,14.503
5001,304383.335,4239210.993,14.274
5002,304387.942,4239207.943,14.070
5003,304386.662,4239208.707,14.649
5004,304384.396,4239210.264,14.658
5005,304382.937,4239215.484,13.127
5006,304391.427,4239199.615,14.625
5007,304384.605,4239215.193,13.790
5008,304390.320,4239210.924,13.144
5009,304390.675,4239211.730,13.527
5010,304384.028,4239214.456,13.600
5011,304389.808,4239213.925,14.544
5012,304386.617,4239215.258,14.364
5013,304391.342,4239213.373,14.222
5014,304387.307,4239215.455,14.641
5015,304392.840,4239212.544,14.665
5016,304393.434,4239213.051,14.460
5017,304398.104,4239211.311,14.500
5018,304392.086,4239231.001,14.012
5019,304400.437,4239226.271,14.478
5020,304394.035,4239229.973,14.158
5021,304398.668,4239227.504,14.207
5022,304395.202,4239229.513,14.610
5023,304397.281,4239228.193,14.747
5024,304413.300,4239239.443,14.534
5025,304393.100,4239237.360,13.732
5026,304406.513,4239240.325,14.140
5027,304397.282,4239236.282,14.241
5028,304404.371,4239241.645,14.766
5029,304398.859,4239236.511,14.614
5030,304402.167,4239242.920,14.585
5031,304401.193,4239243.205,14.177
5032,304400.646,4239251.586,13.631
5033,304418.313,4239250.066,14.883
5034,304422.754,4239257.478,14.703
5035,304444.362,4239229.472,15.309
5036,304441.297,4239242.967,15.053
5037,304438.475,4239245.447,15.158
5038,304414.531,4239252.707,14.844
5039,304422.651,4239256.339,14.656
5040,304403.724,4239260.289,14.472
5041,304407.046,4239259.478,14.485
5042,304396.572,4239266.052,14.095
5043,304392.975,4239268.788,13.868
5044,304384.229,4239274.796,13.486
5045,304381.430,4239277.784,13.361
5046,304383.257,4239280.128,13.392
5047,304385.887,4239282.551,13.315
5048,304383.247,4239285.375,13.283
5049,304399.191,4239274.070,14.077
5050,304399.030,4239273.211,13.864
5051,304412.882,4239264.036,14.686
5052,304410.092,4239266.820,14.294
5053,304409.109,4239270.575,14.065
5054,304404.520,4239284.708,13.264
5055,304413.820,4239280.421,13.877
5056,304405.821,4239281.551,13.849
5057,304412.027,4239270.964,14.330
5058,304411.379,4239292.844,13.609
5059,304419.455,4239287.448,14.028
5060,304411.155,4239297.672,13.137
5061,304416.666,4239275.539,14.303
5062,304419.244,4239271.756,14.711
5063,304417.267,4239273.186,14.666
5064,304429.346,4239271.806,13.916
5065,304395.427,4239196.949,14.720
5066,304385.302,4239202.883,14.375
5067,304376.617,4239207.803,13.768
5068,304387.468,4239182.638,14.644
5069,304383.638,4239175.801,14.600
5070,304365.398,4239189.157,13.613
5071,304376.336,4239188.549,14.457
5072,304374.562,4239188.912,14.751
5073,304372.366,4239190.320,14.687
5074,304370.105,4239191.848,13.683
5075,304364.947,4239173.767,14.812
5076,304362.757,4239175.179,14.624
5077,304362.177,4239175.720,14.239
5078,304374.371,4239162.035,14.878
5079,304350.947,4239173.981,13.540
5080,304345.919,4239165.269,13.810
5081,304363.117,4239145.544,14.709
5082,304353.417,4239166.347,13.852
5083,304356.286,4239164.816,14.745
5084,304358.347,4239163.389,14.828
5085,304365.093,4239170.759,14.443
5086,304359.645,4239162.570,14.584
5087,304339.767,4239168.717,13.836
5088,304335.119,4239160.823,13.925
5089,304347.578,4239141.143,14.646
5090,304346.408,4239142.013,14.716
5091,304346.308,4239142.028,15.587
5092,304322.975,4239142.280,13.812
5093,304344.722,4239143.165,15.585
5094,304341.683,4239133.002,15.290
5095,304341.707,4239132.956,14.747
5096,304340.504,4239133.818,14.716
5097,304340.485,4239133.884,15.259
5098,304332.146,4239149.120,13.782
5099,304334.400,4239147.283,13.852

Output ASCII file ready to convert to raster using ArcGIS ASCII to Raster tool:

ncols 12 
nrows 16
xllcorner 304322.975000
yllcorner 4239132.956000
cellsize 10.000000
nodata_value -9999
-9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000
-9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 13.314204 13.489219 13.777728 -9999.000000 -9999.000000 -9999.000000
-9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 13.676529 14.018830 14.154381 14.028621 -9999.000000 -9999.000000
-9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 13.389017 13.669945 13.963754 14.338332 14.328053 -9999.000000 -9999.000000
-9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 13.457683 13.620355 13.675830 14.382545 14.695792 14.690051 -9999.000000
-9999.000000 -9999.000000 -9999.000000 -9999.000000 13.492555 13.541571 13.660249 13.908317 14.387550 14.699870 14.904953 -9999.000000
-9999.000000 -9999.000000 -9999.000000 -9999.000000 13.464333 13.625458 13.731216 14.280302 14.610317 14.831402 15.052487 -9999.000000
-9999.000000 -9999.000000 -9999.000000 13.522997 13.483254 13.698401 14.079552 14.498293 14.710660 14.923066 15.178369 -9999.000000
-9999.000000 -9999.000000 -9999.000000 13.638777 13.534187 13.797617 14.354722 14.583508 14.863345 15.061120 -9999.000000 -9999.000000
-9999.000000 -9999.000000 -9999.000000 13.561675 13.585120 14.176141 14.570483 14.779832 15.017311 -9999.000000 -9999.000000 -9999.000000
-9999.000000 -9999.000000 13.763785 13.526557 13.771358 14.485951 14.650874 14.973502 -9999.000000 -9999.000000 -9999.000000 -9999.000000
-9999.000000 -9999.000000 13.736118 13.764651 14.389090 14.622974 14.929693 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000
-9999.000000 13.770209 14.186704 14.640372 14.825527 14.885883 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000
-9999.000000 13.892236 14.308787 14.729322 14.785682 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000
-9999.000000 14.098530 14.584644 14.712380 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000
-9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000
Improve this answer

Your Answer

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

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