Rasterize irregular grid geocube python: Loosing points at higher resolution

119 Views Asked by At

I'm trying to convert a dataframe with an irregular grid to a regular rectilinear grid. I'm currently using the VIIRS hotspot data and the code below should allow people to read in the csv in a colab notebook or from a local machine

!pip install geocube
!pip install geopandas
!pip install xarray

import pandas as pd
import xarray as xr
import numpy as np
import geocube
import geopandas as gpd
from geocube.api.core import make_geocube
from geocube.rasterize import rasterize_image
from geocube.rasterize import rasterize_points_griddata

df = pd.read_csv('https://firms.modaps.eosdis.nasa.gov/content/notebooks/sample_viirs_snpp_071223.csv')

I got a subset for the data

# Get subset
df_subset = df[(df['longitude'] >= 110) & (df['latitude'] >= -55) & 
               (df['longitude'] <= 180) & (df['latitude'] <=-10)].copy()

and converted this subset to a geopandas dataframe which allowed to specify the projection

gdf = gpd.GeoDataFrame(df_subset, 
                       geometry=gpd.points_from_xy(df_subset.longitude, 
                                                   df_subset.latitude), 
                       crs='EPSG:4326'
                       )
# Convert the 'acq_date' column to datetime format
gdf['acq_date'] = pd.to_datetime(gdf['acq_date'])

But then how do I go about when I want to convert the point data to a raster? I tried using geocube which looks fine when I choose a coarse resolution

from functools import partial
from geocube.rasterize import rasterize_image, rasterize_points_griddata,MergeAlg

from geocube.api.core import make_geocube

geo_grid = make_geocube(
    vector_data=gdf,
    measurements=['frp'],
    resolution=(-1, 1),
    rasterize_function=partial(rasterize_image, 
                               all_touched=True),
)
import matplotlib.pyplot as plt
plt.imshow(geo_grid.frp,vmin=0,vmax=5)

enter image description here

But then when I increase the resolution I loose a lot of points

geo_grid = make_geocube(
    vector_data=gdf,
    measurements=['frp'],
    resolution=(-0.25, 0.25),
    rasterize_function=partial(rasterize_image, 
                               all_touched=False))

enter image description here

How can I rasterize the points to a 0.25 grid where each grid cell shows the sum of all data points within the gridcell?

0

There are 0 best solutions below