2

Following this post, I'ms trugglig wwith interpolation and projection of set of points. I have used griddata in order to interpolate, and seems like interpolation hapenns, however, the results look tilted ,like transform is wrong, but I couldn't find the mistake in the calculation.

This is the original set of points:

enter image description here

This is how I do the interpolation:

points3d = gpd.read_file('shape/clay.shp')
points3d = gpd.GeoDataFrame(
    points3d, geometry=gpd.points_from_xy(points3d.Longitude, points3d.Latitude))
print('number of points:{}'.format(len(points3d)))
points3d = points3d.set_crs(epsg=4326)

totalPointsArray = np.zeros([points3d.shape[0],3])


for index, point in points3d.iterrows():
    pointArray = np.array([point.geometry.coords.xy[0][0],point.geometry.coords.xy[1][0],point['Data']])
    totalPointsArray[index] = pointArray


x=totalPointsArray[:,0]
y=totalPointsArray[:,1]
z=totalPointsArray[:,2]


xmin = min(x)
xmax = max(x)
ymin = min(y)
ymax = max(y)


number of pixels with 1m resolution


nx = (int(xmax - xmin + bbox_size[0]))
ny = (int(ymax - ymin + bbox_size[1]))


xi = np.linspace(xmin, xmax, nx)
yi = np.linspace(ymin, ymax, ny)
xi, yi = np.meshgrid(xi, yi)


zi = inp.griddata((x, y), z, (xi, yi), method='linear')
plt.imshow(zi)

The result seems to be tilted: enter image description here

I couldn't find the error I have in the calculation. My end goal is to have interpolation aligned to the points and not "tilted".

enter image description here

ReutKeller
  • 2,139
  • 4
  • 30
  • 84

1 Answers1

2

In the end the problem was mismatch between geography (lower left) and graphics (upper left). I have gotten my bounding box coordinates (in this case, I have used sentinel-hub bbox , hence I have used the get_polygon function:

bbox.get_polygon()

and then based on the coordinates of the bounding box I have carefully selected the uper left instead of the lower left coordinates of the y-axis.

ReutKeller
  • 2,139
  • 4
  • 30
  • 84