extract

[Frederik Noack]

Dear Tomislav,

I am trying to relate incomes of 8000 poor rural households in 350 villages in the tropics to soil data. I have the coordinates for the villages in csv format and would like to extract soil properties from the ISRIC 1km soil grid for each location in R. I have used the extract function with netcdf files before but I have difficulties to use it with the soil data. Could you please help me with an approach to extract soil properties for each village location in a 10km radius?

Cheers
Frederik

[Tomislav Hengl]

If you need to obtain values of SoilGrids for e.g. <500 points then the most efficient thing is to use the REST.SoilGrids function (http://gsif.r-forge.r-project.org/REST.SoilGrids.html). First test out the REST service by manually adding the coordinates and variables of interest (e.g. http://rest.soilgrids.org/query?lon=5.39&lat=51.57&attributes=ORCDRC), then use the example from REST.SoilGrids and adjust / extend.

However, if you wish to obtain aggregated values of soil properties within a buffer e.g. 10 km around the points of interest, then I would advise that you:

1. Dowload the GeoTiff of interest (bounding box + 10 km in each direction) using the SoilGrids WCS (see http://gsif.isric.org/doku.php?id=wiki:tutorial_soilgrids#wcs_data_access i.e. this code -> https://code.google.com/p/gsif/source/browse/trunk/SoilGrids1km/tutorial_SoilGrids.R)
   
2. Import the GeoTiffs using the raster package.
3. Buffer the values by setting up the appropriate radius (see e.g. http://rossijeanpierre.wordpress.com/2014/01/31/extracting-circular-buffers-from-a-raster-in-r-12/ but this can also be done via the extract function I think).

Try to make some R code and then send me R code chunks if you experience any problems. If you wish to share files / scripts with me use either DropBox (username: "tom....@gmail.com") or similar. Please do not send data via this mailing list.

HTH and apologies for some delay in my reply,

T. (Tom) Hengl

Researcher @ ISRIC - World Soil Information

Url: http://www.wageningenur.nl/en/Persons/dr.-T-Tom-Hengl.htm

Network: http://profiles.google.com/tom.hengl

Publications: http://scholar.google.com/citations?user=2oYU7S8AAAAJ

[Jorge de Jesus]

Hi to all

Just another suggestion, but in this case you have to  do some extra programming.

The REST interface can also inform  what tile you need to donwload for a specific point e.g::

http://rest.soilgrids.org/stream/tile_with_point?lon=5.39&lat=51.57 

This for example will reply with:

T523

With this Tile ID you can again request for the service to  give you the tile

http://rest.soilgrids.org/stream?tiles=T523

The query can also include parameters specific for parameters, e.g:

http://rest.soilgrids.org/stream?tiles=T523&attributes=CEC,BLD

Just pay attention that you can't request more than 5 tiles at the same time or the service will block you

From the download tiles you run the buffer code.

Personally, I think the solution from Tom is clear and easier to implement.

Jorge

--
You received this message because you are subscribed to the Google Groups "Global Soil Information" group.
To unsubscribe from this group and stop receiving emails from it, send an email to global-soil-infor...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

OpenPGP Key: 0x7212572C

[Frederik Noack]

Dear Tom and Jorge,

Thanks a lot for your answers! I tried

village=data.frame(GEO[,c("gvillcode","Long.dec.deg","Lat.dec.deg")])

names(village)=c("id","lon","lat")

coordinates(village)= ~lon+lat

proj4string(village) = CRS("+proj=longlat +datum=WGS84")

soil=REST.SoilGrids(c("ORCDRC","PHIHOX","SNDPPT","SLTPPT","CLYPPT","GRAVOL","CEC","DBR"))

soil_village=over(soil, village)

 

where GEO is the file with my village coordinates in decimal degrees. I think, I have now matched the village coordinates with the soil data of that specific grid. However, since I believe that the fields of the households for agricultural production are mainly located outside this grid I would like to use a buffer. My villages are in 20 tropical countries. How can I use the proposed code for an area defined by coordinates and not countries?

> wg.url <- url("http://gsif.isric.org/lib/exe/fetch.php?media=admin.af.rda")

> load(wg.url)

> proj4string(admin.af) <- "+proj=longlat +datum=WGS84"

> country.af <- as(admin.af, "SpatialLines")

> ghana <- admin.af[admin.af$FORMAL_EN=="Republic of Ghana",]

 

Cheers

Frederik

--
You received this message because you are subscribed to a topic in the Google Groups "Global Soil Information" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/global-soil-information/HwzLwTbtTv4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to global-soil-infor...@googlegroups.com.