[GIS] Continuous smoothing of raster using R

graphicsrrastersmoothingweighted-overlay

I have a dataset on some villages like this on a grid:

I want to fill the space of the whole grid so that each grid cell has a value that reflects the pressure from the surrounding villages (i.e., the value would be something like average distance to all villages, weighted by the size of each village).

As a first attempt, I tried to smooth it with:

    villageRaster <- focal(villageRaster, w=matrix(1, 9, 9), mean)

to give:

But there are still large white areas with no values.

Can someone suggest a nice way to smooth across the whole space?

I am using R for this preferably.

Best Answer

You could look into this approach as suggested Smoothing raster map using R?:

Here are some ideas.

With base plot you can do

 plot(x, interpolate=TRUE)

You can also resample your data using disaggregate (x, 5, method='bilinear')

 y <- disaggregate x, 5, method='bilinear')

Or indeed smooth it using a focal operation

 y <- focal(x, w=matrix(1, 5, 5), mean)

Or a combination

 y <- disaggregate(x, 5)
 y <- focal(y, w=matrix(1, 5, 5), mean)

The question whether doing this is a good idea or not is another matter, that I'll leave to you to decide

Related Solutions

[GIS] How to iterate through every cell in a continuous raster

I see this has already been solved by the Original Poster (OP), but I'll post a simple solution in python just in case anyone in the future is interested in different ways to solve this problem. I'm partial to open source software, so here's a solution using GDAL in python:

import gdal

#Set GeoTiff driver
driver = gdal.GetDriverByName("GTiff")
driver.Register()

#Open raster and read number of rows, columns, bands
dataset = gdal.Open(filepath)
cols = dataset.RasterXSize
rows = dataset.RasterYSize
allBands = dataset.RasterCount
band = dataset.GetRasterBand(1)

#Get array of raster cell values.  The two zeros tell the 
#iterator which cell to start on and the 'cols' and 'rows' 
#tell the iterator to iterate through all columns and all rows.
def get_raster_cells(band,cols,rows):
    return band.ReadAsArray(0,0,cols,rows)

Implement the function like this:

#Bind array to a variable
rasterData = get_raster_cells(band,cols,rows)

#The array will look something like this if you print it
print rasterData
> [[ 1, 2, 3 ],
   [ 4, 5, 6 ],
   [ 7, 8, 9 ]]

Then, iterate through your data with a nested loop:

for row in rasterData:
    for val in row:
        print val
> 1
  2
  3
  4...

Or maybe you want to flatten your 2-D array with a list comprehension:

flat = [val for row in rasterData for val in row]

Anyways, while iterating through the data on a cell-by-cell basis its possible to throw some conditionals into your loop to change/edit values. See this script I wrote for different ways to access the data: https://github.com/azgs/hazards-viewer/blob/master/python/zonal_stats.py.

LiDAR to DEM – How to Determine Average Points Spacing and Cell Size in ArcGIS

It really depends on what you're going for here. If your goal is to make accurate measurements, you will want the grid spacing to be smaller, around the average point spacing of the dataset. If this is part of a process of creating orthophotos or analyzing watersheds, or any other task which requires a more generalized terrain model, you will want the spacing to be larger.

Also, if you've filtered the data to the ground points, USE THOSE! You can ignore the other classifications.

Three times the size of the point spacing is generally considered a good idea for generalization because it allows the grid cell three points to average out to a theoretically consistent value. In my profession, sometimes when I produce orthophotos I allow a huge grid size, often 8-12 times the size of the point spacing.

My point (heh) is: there is no right or wrong way to produce a bare earth model. Sure, there are specifications sometimes (and if this is one of those cases, I suggest you disregard whatever I'm saying and follow the specifications), but producing what you need to suit your problem can be totally heuristic.

Do you want a more rigid, presumably more accurate model with harder edges on which you can model terrain more precisely? Use a smaller spacing. Do you want a more general terrain model that is smoother, less precise, but shows the overall terrain characteristics of a large area? Use a large spacing.

A lot of this is trial and error. That's the beauty of this kind of thing: it's about 30% science and 70% art.

Best Answer

Related Solutions

[GIS] How to iterate through every cell in a continuous raster

LiDAR to DEM – How to Determine Average Points Spacing and Cell Size in ArcGIS

Related Question