0

I am new in GIS, and I downloaded a gridded dataset from this source

The data come with a raster file as well as a range of attribute data files. In ArcGIS I can easily join the .mdb Microsoft Access attribute table with the .adf raster. However, I am now in a Mac, that can't handle .mdb. Fortunately the data is also downloadable in an interchangeable format, including a .tiff raster file together with several .txt files with attributes.

How can I link the TIFF with the TXT data in R or QGIS from a Mac?

Once uploaded, how can I select one particular attribute (and discard the others) similar to what I would do with the lookup function? The idea is to finally save a raster in tiff with only one specific column among all the different variables included in the attribute tables.

  • What are the rows and columns of the attribute table, and how do they relate to the values in the grid? Can you show a sample of the attribute table? Are there values in the grid correspond to columns in the attribute table? – Spacedman Jul 5 '17 at 22:47
  • (Your linked sample data is a bit big really...) – Spacedman Jul 5 '17 at 22:52
  • From the ReadMe file: The zip file includes a raster GIS file (WISE30sec) as well as a range of soil attribute data files; Using ArcMap®, users may join the raster data to the derived soil properties files, as managed the MSAccess® database (WISE30sec.mdb). Linkage is through the map unit code or grid cell identifier (NEWSUID) of the raster set and the NEWSUID of the various soil attribute data files. The above files are also provided in TIFF and TXT format (see folder “...\Interchangeable_format”) so that they may be easily accessed using other systems. – fede_luppi Jul 6 '17 at 6:31
  • The columns are different soil properties: pH, nitrogen content, C:N ratio, etc. I guess the rows are polygons, – fede_luppi Jul 6 '17 at 6:33
2

The data appears to be in three formats: a .mdb in ./Data/, a .tif and some txt files in ./Interchangeable_format/, and an Arc/Info binary coverage in ./GISfiles with a raster layer called wise30sec_fin.

Both the tif and the Arc/Info coverage can be read by R or QGIS via the underlying GDAL library. The tif appears to be just integers that will need to be used for lookups, but the Arc/Info coverage includes those and is a bit richer too, so I'll use that...

Use the raster package to read it in:

> library(raster)
> r = raster("./wise30sec_fin/")
> r
class       : RasterLayer 
dimensions  : 16753, 43201, 723746353  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -180, 180.0083, -55.98333, 83.625  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /home/rowlings/Downloads/wise_30sec_v1/GISfiles/wise30sec_fin 
names       : wise30sec_fin 
values      : 0, 32156  (min, max)
attributes  :
          ID MU_GLOBAL MUGLB_NEW COVERAGE CODE CLIMATE_CO    NEWSUID     COUNT
 from:     0         0         0        0             Y-0            502784485
 to  : 32156     31111     32156        5    C    C-32156 WD30032156        51
 VALUE30SEC
          0
      32156

Note this hasn't actually read it into memory because of its size. R will get bits of it as needed, and we have to be careful not to read it all into memory unless we have a lot of memory. You may struggle to process the whole raster in R, so crop it to your study area if you have one.

Note the raster is a single band, and has a set of attributes. This is a data frame that links integer values in the raster to attribute values. Lets look at the first few, by using levels on the raster. The raster is a single band, so we get the first attribute table by extracting with [[1]]:

> levels(r)[[1]][1:4,1:5]
  ID MU_GLOBAL MUGLB_NEW COVERAGE CODE
1  0         0         0        0     
2  2         2         2        4    A
3 36        36        36        4    B
4 37        37        37        4    A

So to get the information for a location, get the value of the raster at that location and that's the ID in the attributes table:

> r[5000,7000]

4713 
> subset(levels(r)[[1]], ID==4713)
       ID MU_GLOBAL MUGLB_NEW COVERAGE CODE CLIMATE_CO    NEWSUID COUNT
1999 4713      4713      4713        4    D     D-4713 WD40004713  8522
     VALUE30SEC
1999       4713

Now I can read in one of the data files and look up in that with the NEWSUID code:

> data = read.table("./Interchangeable_format/HW30s_wD1.txt",head=TRUE,sep=",")

The data is one row shorter then the attribute table, so its possible one NEWSUID code is missing:

> dim(data)
[1] 16413    55
> dim(levels(r)[[1]])
[1] 16414     9

and its the empty string which I guess is code for the sea.

Anyway, we can now select the values from the data file since we know the NEWSUID value at r[5000,7000]:

> subset(data, NEWSUID=="WD40004713")
         NEWSUID NofComponents SoilUnits Layer PROP_aw MiscUnits PROP_misc
12604 WD40004713             1   GLm100     D1     100                   0
      Drain DrainProp DrainNum DrainMin DrainMax TopDep BotDep CFRAG CFRAG_std
12604     P       100        2        2        2      0     20     1         1

So that's the fundamentals of looking up data values given the raster and some data. You can do all sorts of clever things like adding the data to the levels(r) but I think the best thing to do might depend on what you actually want to end up doing with the data - if its just sampling at a few points you might want to do something different to if you want to create a whole load of global maps of a quantity in the data files.

There's a useful answer on manipulating raster level attributes here:

https://stackoverflow.com/questions/28617678/using-a-raster-attribute-from-a-multi-attribute-raster-for-colour-levels-in-a-pl

but watch out - some of those things might try and load the whole raster into memory and that might be pretty slow or hang your computer.

Alternatively, using exclusively the files in ./Interchange_files, the tif file stores pixel values, which are in the .tsv file, where you can lookup NEWSUID values to look up in the data file. Again, the large size stops me wanting to do this for the whole raster but everything is well-named.

  • Dear Spacedman, I have now been able to look at your response. First of all, thanks for your help. The final output I need is a global map of soil C:N ratio. In this case, I would want to somehow join the raster with only one particular column from data, "CNrt", plot it, and save this C:N map as a tiff. How could I do this? I have tried new <- subs(r, data, by="NEWSUID, which="CNrt") but it seems to take ages to complete, and I don't even know if this is the way to go, so I had to cancel. Please, how can I get a map of CNrt only? Thanks – fede_luppi Jul 26 '17 at 8:46
  • I finally let it run. After a while the process finished with an empty layer – fede_luppi Jul 26 '17 at 15:25
  • please, can you indicate if you need further information to be able to help me? – fede_luppi Aug 1 '17 at 9:52

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.