0

I have one raster that is categorical values for land cover classes (from the NLCD) at 30m resolution. I also have another raster for land surface temperature (LST) data at 70m resolution. I would like to find the dominant NLCD class in each LST pixel. For example, if there are 4 NLCD pixels in a single LST pixel, I would like to find the NLCD class with the most occurences.

I would prefer Python, but would accept an answer in R (since it has become very popular for Raster analysis).

Test dataset located at: https://github.com/arojas314/data-sharing/blob/main/nlcd2019_nyc.zip

2 Answers 2

0

Here is an approach with R

library(terra)
nlcd = rast("nlcd2019_nyc.tif")
lst = rast("lst_scaled_temp.tif")

# create a template that matches the extent of lst and 
# approximately the spatial resolution of nlcd
x = disagg(rast(lst), 2)

# nearest neighbor resample nlcd to the template
rnlcd = resample(nlcd, x, "near")

# aggregate to get the dominant value of each lst cell
arnlcd = aggregate(rnlcd, 2, "modal")
arnlcd
#class       : SpatRaster 
#dimensions  : 851, 983, 1  (nrow, ncol, nlyr)
#resolution  : 0.0006308342, 0.0006308461  (x, y)
#extent      : -74.30966, -73.68955, 40.41139, 40.94824  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#source      : memory 
#color table : 1 
#name        : Layer_1 
#min value   :       0 
#max value   :      95 

You do not provide any background info. I wonder if your goal is to summarize lst values by nlcd class. In that case I think it would make more sense to do:

rlst = resample(lst, nlcd)
z2 = zonal(rlst, nlcd, mean)
head(z2)
#  Layer_1 lst_scaled_temp
#1       0        283.6699
#2      11        289.6679
#3      21        300.0840
#4      22        301.9470
#5      23        304.7881
#6      24        308.4624
1
  • I tend to avoid R since I usually integrate python code into larger programs, but wow this method is easy. Thanks!
    – Clouseau
    Commented Jun 2, 2022 at 1:15
0

you have to co-register two data. since you want dominant value in each cell, i think mode resampling method is suitable for your task. here the python function for co-registering using rasterio.
reference: co-register (code not mine - i just changed resampling method to mode):

from rasterio.warp import reproject, Resampling, calculate_default_transform
import rasterio
def reproj_match(infile, match, outfile):
    """Reproject a file to match the shape and projection of existing raster. 
    
    Parameters
    ----------
    infile : (string) path to input file to reproject
    match : (string) path to raster with desired shape and projection 
    outfile : (string) path to output file tif
    """
    # open input
    with rasterio.open(infile) as src:
        src_transform = src.transform
        
        # open input to match
        with rasterio.open(match) as match:
            dst_crs = match.crs
            
            # calculate the output transform matrix
            dst_transform, dst_width, dst_height = calculate_default_transform(
                src.crs,     # input CRS
                dst_crs,     # output CRS
                match.width,   # input width
                match.height,  # input height 
                *match.bounds,  # unpacks input outer boundaries (left, bottom, right, top)
            )

        # set properties for output
        dst_kwargs = src.meta.copy()
        dst_kwargs.update({"crs": dst_crs,
                           "transform": dst_transform,
                           "width": dst_width,
                           "height": dst_height,
                           "nodata": 0})
        print("Coregistered to shape:", dst_height,dst_width,'\n Affine',dst_transform)
        # open output
        with rasterio.open(outfile, "w", **dst_kwargs) as dst:
            # iterate through bands and write using reproject function
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=dst_transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.mode) # mode resampling method

I've checked and it seems working.

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.