Creating raster file from point data

Coverting point data shapefile to a raster file.¶

In [1]:

import geopandas as gpd

gdf=gpd.read_file("Synthetic_Subsidence_Data.shp")

This code imports the Python package geopandas and reads in a shapefile named Synthetic_Subsidence_Data.shp using the read_file function from geopandas.

What is shapefile? A shapefile is a popular geospatial vector data format used to store geographic data. It contains information such as geographical features, coordinates, and attributes of the data.

By using geopandas, the shapefile data can be read as a GeoDataFrame, which is a pandas DataFrame that includes a column of geospatial data, usually represented as a series of points, lines or polygons.

Therefore, the resulting df object is a GeoDataFrame that contains the data and geometries (points, lines or polygons) from the Synthetic_Subsidence_Data.shp shapefile.

To see this tabular data you can use: head()

In [2]:

gdf.head()

Out[2]:

	Unnamed_ 0	Latitude	Longitude	Subsidence	geometry
0	0	-9.995486	64.556974	0.001131	POINT (64.55697 -9.99549)
1	1	-16.447512	77.177932	0.000139	POINT (77.17793 -16.44751)
2	2	8.360813	59.225169	0.002132	POINT (59.22517 8.36081)
3	3	-6.795347	106.702294	0.000019	POINT (106.70229 -6.79535)
4	4	12.859304	40.289230	0.000143	POINT (40.28923 12.85930)

Visualising Point Data¶

To learn more about visualisation shapefile you can follow this tutorial

In [3]:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt


# Create a figure and axis object
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection=ccrs.PlateCarree()))

# Add coastlines and features to the plot
ax.coastlines(resolution='10m')
ax.scatter(gdf['Longitude'],gdf['Latitude'],c=gdf['Subsidence'])
# Set the extent of the plot to the Indian Ocean region
ax.set_extent([40, 110, -20, 30], crs=ccrs.PlateCarree())

# Show the plot
plt.show()

This code is used to create a map plot using the Python packages cartopy, matplotlib, and geopandas. It visualizes data about subsidence in the Indian Ocean region.

Here is a breakdown of what the code does:

First, the necessary packages are imported, including cartopy.crs for coordinate reference systems, cartopy.feature for map features, and matplotlib.pyplot for plotting.
Then, a figure and axis object are created using the plt.subplots() function. The figsize parameter is used to set the size of the figure in inches, and the projection parameter sets the coordinate reference system to Plate Carrée (longitude and latitude).
Next, the ax.coastlines() function adds coastline features to the plot. The resolution parameter is used to set the detail level of the coastlines.
The ax.scatter() function adds data points to the plot. It takes the longitude and latitude columns from the gdf GeoDataFrame as inputs, and uses the c parameter to color the points based on the Subsidence column.
The ax.set_extent() function sets the extent of the plot to the Indian Ocean region, which is defined by the range of longitude and latitude values provided in the list [40, 110, -20, 30]. The crs parameter is used to specify the Plate Carrée coordinate reference system.
Finally, the plot is displayed using plt.show()

Interpolating Point Data to Raster¶

This code performs interpolation of point data onto a regular grid and writes the resulting raster to a file using the Python packages numpy, scipy, and rasterio

Extracting Relevant Data¶

The first section of the code extracts the relevant data from the gdf GeoDataFrame, including the subsidence values, longitude, and latitude.

Creating the Grid¶

The next section of the code sets the resolution of the output raster and creates a regular grid of points using numpy.meshgrid().

Interpolating Data¶

The scipy.interpolate.griddata() function is used to interpolate the point data onto the regular grid. The method parameter is set to ‘cubic’ to use cubic interpolation.

Defining the Transformation and CRS¶

The Affine class from rasterio.transform is used to define the transformation matrix for the output raster.

Writing file as TIF¶

The interpolated data is then written to a GeoTIFF file using rasterio. First, an Affine transformation and CRS (coordinate reference system) are defined based on the grid specifications. The interpolated data is then written to the file using the rasterio.open() function, specifying the file name, dimensions, data type, CRS, and transformation information. The write() method is used to write the interpolated data to the file, and the file is closed using the close() method.

In [4]:

import numpy as np
from scipy.interpolate import griddata
import rasterio
from rasterio.transform import Affine
from rasterio.crs import CRS

# Extract relevant data
subsidence = gdf['Subsidence'].values
lon = gdf['Longitude'].values
lat = gdf['Latitude'].values

# Set resolution and create grid
resolution = 0.1
x_range = np.arange(lon.min(), lon.max()+resolution, resolution)
y_range = np.arange(lat.min(), lat.max()+resolution, resolution)
grid_x, grid_y = np.meshgrid(x_range, y_range)

# Interpolate data onto grid
points = list(zip(lon, lat))
grid_z = griddata(points, subsidence, (grid_x, grid_y), method='cubic')

# Set negative values to 0
grid_z[grid_z < 0] = 0

# Define transformation and CRS
transform = Affine.translation(grid_x[0][0]-resolution/2, grid_y[0][0]-resolution/2) * Affine.scale(resolution, resolution)
crs = CRS.from_epsg(4326)

# Write interpolated raster to file
interp_raster = rasterio.open('Subsidence.tif',
                              'w',
                              driver='GTiff',
                              height=grid_z.shape[0],
                              width=grid_z.shape[1],
                              count=1,
                              dtype=grid_z.dtype,
                              crs=crs,
                              transform=transform)
interp_raster.write(grid_z, 1)
interp_raster.close()

Plotting raster file¶

To learn more about plotting a raster you read this tutorial.

In [5]:

import rioxarray

subsidence_raster=rioxarray.open_rasterio("Subsidence.tif")

fig=plt.figure(figsize=[12,8])
ax = fig.add_axes([0,0,1,1],projection=ccrs.PlateCarree())
raster_image=subsidence_raster.plot(ax=ax,cmap='jet')
ax.coastlines(resolution='10m')
plt.show()