Creating a lat/lon grid from lat/lon point data of a given variable in python?

user_18568199 0
Viewed 668

I have a large dataframe, df, which I have created from multiple files of scattered, irregular data.

The df is very long (8131596) and has the columns: date, latitude, longitude, var1, var2.

I would ideally like to create a 2D grid/map of certain latitude/longitude bounding box, and apply my df so any data that falls within that lat/lon is included. Or, create 2D lat/lon grid from my data and I can cut it later (ie, end up with a 2D 'image' / matrix).

I will then interpolate between these irregular points, ending up with essentially a data field for each of the variables in the df.

However I am a bit stuck with what is the best way to go about this. I think perhaps mesh grid will be involved, but I can't seem to find any previous questions/help that relates this situation.

Any suggestions greatly appreciated!

pandas geospatial coordinates interpolation python
1 Answers

This is for var1 you would have to repeat this for var2

import numpy as np
import pandas as pd
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import rasterio
from rasterio.crs import CRS

# define interpolation inputs
points = list(zip(df.longitude,df.latitude))
values = df.var1.values # change to var2 for second grid

# define raster resolution
rRes = 50

# create coord ranges over the desired raster extension
xRange = np.arange(df.longitude.min(),df.longitude.max()+rRes,rRes)
yRange = np.arange(df.latitude.min(),df.latitude.max()+rRes,rRes)

# create arrays of x,y over the raster extension
gridX,gridY = np.meshgrid(xRange, yRange)

# interpolate over the grid
gridPh = griddata(points, values, (gridX,gridY), method='linear')

# show interpolated values
plt.imshow(gridPh)
plt.colorbar()

# definition of the raster transform array
from rasterio.transform import Affine
transform = Affine.translation(gridX[0][0]-rRes/2, gridY[0][0]-rRes/2)*Affine.scale(rRes,rRes)
transform

# get crs as wkt
# use your crs here
rasterCrs = CRS.from_epsg('32718')
rasterCrs.data

#definition, register and close of interpolated raster
interpRaster = rasterio.open('./interpRaster3.tif',
                                'w',
                                driver='GTiff',
                                height=gridPh.shape[0],
                                width=gridPh.shape[1],
                                count=1,
                                dtype=gridPh.dtype,
                                crs=rasterCrs,
                                transform=transform,
                                )
interpRaster.write(gridPh,1)
interpRaster.close()
user_127754320
Related