I have percent sand/silt/clay values for 30 sample points on grid. I used these values to create 3 interpolated IDW rasters. I want to somehow (I think in raster calculator) to make a code using if/then statements to take values from each IDW and output a texture.

Example: if raster "sand_IDW" value > 50 and < 70 and raster "clay_IDW" value > 20 and <40 then, "Texture" = Sandy clay.

  • 1
    Unfortunately it's probably not going to be that easy. You may want to take a look at gis.stackexchange.com/questions/141415 (which I'm pretty sure this is a duplicate of though you've done the interpolation) and gis.stackexchange.com/questions/24501 as well as their linked questions. And for the record you should probably specify what software you're using. – Chris W Nov 17 '15 at 20:11
  • If you're familiar with the OGR/GDAL library, you could probably use the ogr2ogr command line tool with the statement -where "Sand_IDW between 50 AND 70. See gdal.org/ogr2ogr.html and w3schools.com/sql/sql_where.asp – geoeye Nov 17 '15 at 23:04
  • The simplest way, assuming you are using ArcGIS is to compute all 4 rasters first and find cell statistics. 1st raster Con(sand_IDW>50,1,0), 2nd raster Con(sand_IDW<70.1,0)... Areas where cell statistics (SUM)=4 is Sandy Clay. etc – FelixIP Nov 18 '15 at 0:29
  • @FelixIP Pretty sure that method wouldn't work. Multiple conditions would have the same sum. Consider the three variables sand clay loam. Values of 20/40/40, 40/20/40, and 40/40/20 all sum to 100, yet they'd be very different textures. The next step is maintaining identity with the value, which is what the visible answer at my first link does. Then there's the IF tree method this asker seeks, which is also there but deleted so you can't see it without sufficient rep. – Chris W Nov 20 '15 at 18:13
  • @ChrisW sand, clay =60,30. Con (sand>50,1,0)+con (sand <70,1,0)+con (clay>20,1,0)+con (clay <40,1,0) = 4. I am suggesting totals of con rasters not original – FelixIP Nov 20 '15 at 18:38

Grass 7 has an add-on called r.soils.texture that does compute soil texture from a %clay and %sand layer.

However, it did not work on my computer (yet).

| improve this answer | |
  • 1
    Can you be more specific on what didn't work? Was it that the add-on didn't run at all, the output was incorrect, etc. – papadoo Nov 30 '15 at 17:18

Dunno why this got bumped, but it got me thinking, so.

If you're going to classify in Arc, you probably need (at minimum) a lookup table with columns clay|silt|sand|texture. Each row would need to hold every possible combination of three positive integers that add to 100%, and then a class for that combo, defined with reference to an existing texture triangle. Then you'd have to round all your rasters off to integer before using the lookup table in raster calculator. Higher precision would require a much longer table, not that you'd bother going past 2 decimal places. Its still a limited, blunt-force approach.

I have a workflow for a quick'n'dirty texture map using R on my blog, if you're happy to step away from ArcGIS. It relies on a package that has been written for exactly this purpose, the underpinnings of which are far more sophisticated than a lookup table (mad props to the author, Julien Moeys).

Your input rasters will need to be correctly projected and in a format that GDAL can handle.

Short version:


clay_src <- 'path/to/clay/raster'
silt_src <- 'path/to/silt/raster'
sand_src <- 'path/to/sand/raster'
inputs <- c(clay_src, silt_src, sand_src)

# read in rasters, stack and promote to SpatialPixelsDataFrame
SSC <- stack()
for (i in 1:length(inputs)) {
  rn  <- inputs[i]  
  r   <- raster(rn)
  SSC <- addLayer(SSC, r)

SSCP <- as(SSC, 'SpatialPixelsDataFrame')

The SpatialPixelsDataFrame is cool because it lets you hang an attribute list off each pixel and treats those lists like a table. Everything stays organised and in sequence as you do your processing. Downside is that it all sits in memory so there are some size limitations.

# do some tidying up

names(SSCP@data) <- c('CLAY', 'SILT', 'SAND')
SSCP@data <- round(SSCP@data, 2)
SSCP@data$raw_totals <- rowSums(SSCP@data[, 1:3])

# this normalises the three datasets so they stay proportional but add to 100%:

SSCP_norm <- TT.normalise.sum(tri.data = SSCP@data, residuals = T)

definitely check the residuals. If you find that too much of your dataset is outside 95-105% you've got a problem with your kriged rasters - the bigger the variation, the more of a problem it is. If they look ok:

# append the normalised data to the original set - its good practice to keep both

colnames(SSCP_norm)[1:3] <- paste0(colnames(SSCP_norm)[1:3], "_n")
SSCP@data <- cbind(SSCP@data, round(SSCP_norm, 2)) 

# the following produces a new column with the classification for each pixel (pretty much!)
# check the soiltexture vignette to pick out a texture triangle suited to your region 
# and dataset (class size limits matter!), there's a bunch pre-defined
SSCP@data <- cbind(SSCP@data, 
                  "TEXCLASS" = TT.points.in.classes(tri.data  = SSCP@data[, c('CLAY_n', 'SILT_n',  'SAND_n')],
                                                    css.names = c('CLAY_n', 'SILT_n', 'SAND_n'),
                                                    class.sys = "AU2.TT", 
                                                    PiC.type  = "t",
                                                    collapse  = ', ')

the output may need some tidying up to handle data that fell right on the edge of >1 texture class, but after that its pretty easy to convert to numeric and export a final classified raster.

| improve this answer | |

Please, try r.soils.texture in grass. inside GRASS 7.0 you can install r.soils.texture with this command line: g.extension extension=r.soils.texture operation=add

you can run the addon r.soils.texture for obtain texture raster file with this command line: r.soils.texture sand=name clay=name scheme=name output=name

You can retrieve example data or schema texture files from this site: http://maplab.alwaysdata.net/soilstools.html

If there are problems please ask me for help.

| improve this answer | |
  • Welcome to GIS SE! As a new user please take the Tour. This doesn't answer the question. Please edit your answer to expand on how to use what you are suggesting, not just a one line "try this". – Midavalo Oct 7 '16 at 19:17
  • I've explaned how to install, the previus user said: "However, it did not work on my computer (yet)." – Gianluca Massei May 1 at 15:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.