0

I have two raster files with the same coordinator. I want to first calculate the difference between the two raster files. Then, I want to calculate the percentage of area where the number increases by 10 or more (the number decreases by 10 or more). I am wondering what is the fastest way to do so. Do I need to polygonize the raster to calculate the percentage or I can do it faster as raster format?

Vince
  • 20,017
  • 15
  • 45
  • 64
frituodece
  • 9
  • 1
  • 2
    Stick with raster calculations if possible for the sake of efficiency. Sounds like you could use some simple band math in a raster calculator. For example Raster1 - Raster2 = Difference Raster. Then classify the difference raster so that value >10 or <10 are assigned a value such as 1 and everything else a value of 0. More details: https://gis.stackexchange.com/a/110177/8104 – Aaron Aug 23 '22 at 17:31

1 Answers1

1

Here is a programmatic solution using gdal and pandas.

First using rasterio, numpy, and geocube we can create some sample data :

from osgeo import gdal, ogr, osr
import numpy
import rasterio as rio
from geocube.api.core import make_geocube

this is an existing vector consisting of a single polygon geometry


poly_file = ogr.Open('PDC.shp')


output filename for new raster


output_file = "out_raster.tif"


make geocube from polygon


out_grid = make_geocube(
    vector_data='PDC.shp',
    measurements=["id"],
    resolution=(-25, 25),
    fill=-9999,
)


write original


out_grid['id'].rio.to_raster('in_raster.tif')


load in original raster with gdal


raster_original = gdal.Open('in_raster.tif')


create new raster


driver = gdal.GetDriverByName("GTiff")
raster_out = driver.Create(output_file,
                           len(out_grid.x),
                           len(out_grid.y),
                           1,
                           gdal.GDT_Int16)


create data for new raster and write to output


raster_values = numpy.random.choice(
    [-10, 1, 11], size=(len(out_grid.y), len(out_grid.x)), p=[0.33, 0.34, 0.33])
raster_out.GetRasterBand(1).WriteArray(raster_values)


georeference new raster to match original


srs = osr.SpatialReference()
srs.ImportFromEPSG(3857)
dest_wkt = srs.ExportToWkt()
raster_out.SetProjection(dest_wkt)
raster_out.SetGeoTransform(raster_original.GetGeoTransform())


clear variables to close files


raster_original = None
raster_out = None

Now compute the differences like @Aaron recommended in the comment, categorize them based on the values, then calculate the proportions:

from osgeo import gdal
import pandas as pd

#read in rasters
r1 = gdal.Open('in_raster.tif')
r2 = gdal.Open('out_raster.tif')


#convert to arrays
arr1 = r1.GetRasterBand(1).ReadAsArray()
arr2 = r2.GetRasterBand(1).ReadAsArray()


#compute differences 
out_arr = arr2 - arr1


#convert to categorical and rename categories 
buckets = pd.Categorical(pd.cut(out_arr.flatten(
), [-12, -1, 1, 12])).rename_categories(['fell', 'remained', 'rose'])


#calculate proportions
proportions = buckets.value_counts() / len(buckets)


#print out result
print(proportions)

Output:

enter image description here

Kartograaf
  • 2,902
  • 7
  • 23