2

I am trying to reclassify the CORINE map for Hungary to 4 classes in R. but after reclassification, when I import the asc file to QGIS it loads it empty. Am I doing something wrong?

x<-scan("C:/Users/HP/Documents/Corine_Hungary/Corine_Hungary_2006.asc", sep=" ")
row=length(x)
l <- 1:row
for (i in seq(along=l)){
      if (is.na(x[i])){
      x[i]=0
    } else if (x[i]==0){ 
      x[i] = 0
    }else if (x[i]<4){
      x[i]=1
    }else if (x[i]==5){
      x[i]=1
    }else if (x[i]==4){
      x[i]=2
    }else if (x[i]<10){
      x[i]=2
    }else if (x[i]<30){
      x[i]=3
    }else if (x[i]<35){
      x[i]=2
    }else if (x[i]<45){
      x[i]=4
    }else if (x[i]>45){
      x[i]= "NODATA"
    } 
  }

x=t(x)
write.table(x, "c://Users/HP/Documents/Corine_Hungary
  • 1
    Is your code is finished write.table(x, "c://Users/HP/Documents/Corine_Hungary... at this point? – huckfinn Feb 10 '16 at 9:12
  • write.table(x, "c://Users/HP/Documents/Corine_Hungary/Corine_Reclassified_2006.asc", sep=" ", row.name=FALSE, col.name=FALSE) sorry a result of bad copy paste \ – Mariam Malka Karagulyan Feb 12 '16 at 7:48
  • If the asc file is just an ascii file of numbers you should be able to test your method on a tiny file. Type some numbers into a file, process it as above, see what comes out. – Spacedman Feb 12 '16 at 23:10
  • @spacedman ..ASCII Grid - The ArcInfo ASCII Grid format is an ArcInfo Grid exchange file with a single file extension *.asc ..also possible – huckfinn Feb 12 '16 at 23:33
  • If its a valid ArcInfo ASCII grid then your scan would fail because those files have headers and scan would complain because by default it expects only numbers unless you use the what argument to read characters. And then your code would try and compare the header text with numbers and complain. Have you chopped the header off the file? Open the file in a text editor and say what you see. The first few lines should specify "NCOLS" and "NROWS" and so on. Your code will generate error messages on scan with a valid ArcInfo ASCII Grid file. – Spacedman Feb 13 '16 at 8:56
3

I cannot recover your error, because do not know the content of the Corine_Hungary_2006.asc file. You could either work with the vector stuff you will find under http://www.eea.europa.eu/data-and-maps/data/clc-2006-vector-data-version-3 and select/ classify the right vectors together in a layer. Or you follow this how to:

At first, get the corine data form the EEA website

a. The hungarian grid data from http://www.eea.europa.eu/data-and-maps/data/eea-reference-grids/zipped-shape-files-hungary/zipped-shape-files-hungary/at_download/file unpack the hu.zip into the folder grid-hungary.d.

b. The corine data 2006 from http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-2/clc-2006-100m/g100_06.zip/at_download/file and unpack the g100_06.zip it into the folder grid-europe.d

Cut out the hungarian part out of the very large (59000 x 46000 cells) european grid using the borders of hungarian grid. To get the bounding box I use commands from shapelib http://shapelib.maptools.org/ but a pointer in a gis application will do the job too.

# Read the bounding box
$ shpdump grid-hungary.d/HU_1K.shp  | head
Shapefile Type: Polygon   # of Shapes: 94272
File Bounds: (4786000,2549000,0,0)
         to  (5280000,2896000,0,0)
...

At second, put th things together in a R-Script and do the cut, re-classification and storage jobs. Here is the script using the packages rgdal and gdalUtils:

# load the packages
require("rgdal")
require("gdalUtils")

# Clean up every thing
rm(list = ls())

# setup the datat directory
setwd("~/dev.d/gis-exchange.d/data.d/corine.d");

# South-West corner
sw <- c(4786000,2549000);

# North-east corner
ne <- c(5280000,2896000);

# reassemble the cut window for gdal_translate
# to North-West and South-East corners
cut.win   <- c(sw[1],ne[2],ne[1],sw[2])

rst <- gdal_translate(
               src_dataset='grid-europe.d/g100_06.tif',
               dst_dataset='grid-hungary.d/g100_06.tif',
               projwin=cut.win, output_Raster=TRUE);

# Only some steps of reclassification ..you have to do the rest
values(rst)[is.na(values(rst))] <- 255
values(rst)[values(rst) >   0 & values(rst) <   5]  <- 1
values(rst)[values(rst) >   4 & values(rst) <  10]  <- 2
values(rst)[values(rst) >= 10 & values(rst) < 255] <- 255
# And so on....
# values(rst)[values(rst) > XX & values(rst) < YY] <- ZZ 

# Store the results
writeRaster(rst, filename="grid-hungary.d/g100_06_class.tif", format="GTiff", overwrite=TRUE)

# EOF

Some screenshots for the single steps:

  1. The european corine grid:

CLC 2006 Europe in QGIS

  1. The hungarian part of the european corine grid

CLC 2006 Hungary in QGIS

  1. The "reclassified stuff"

CLC 2006 Hungary reclassified in QGIS

As you see, I've not fully re-classified the dataset, ..your turn now ..enjoy!

  • I think you can do the reclassify in one step with the reclassify function from the raster package. Just build a 3-column matrix of low/high/replacement values and feed it. – Spacedman Feb 13 '16 at 9:03
  • @Spacedman Interesting, can you post some code to enhence the script with the one step reclassification. – huckfinn Feb 13 '16 at 9:13
  • Added as a new answer. The OPs reclassification scheme isn't clear - what happens to 45? Maybe that's not a class in the data. Anyway, if you want to add reclassify to your answer with the right scheme for the data then please do, I'll del my answer. I just can't handle a 70Mb data download right now to try it on CORINE data! – Spacedman Feb 13 '16 at 14:36
  • @Spacedman Thank you, and please deletions. I think it's cool to learn from all of these code snippets, because design, meaning and usage of programs/ code is sometimes a difficult to handle "help" context. – huckfinn Feb 24 '16 at 19:20
2

If you read your data into a raster package object you can use reclassify. Sample data, a grid of 1 to 25:

> x=raster(ncol=5,nrow=5)
> x[]=1:25
> as.matrix(x)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15
[4,]   16   17   18   19   20
[5,]   21   22   23   24   25

Then set up a two-column mapping of old-new values (you can also make a three-column mapping if you want to map within ranges):

> map = cbind(1:25, c(10,10,20,10,10,rep(20,10),rep(30,10)))
> head(map)
     [,1] [,2]
[1,]    1   10
[2,]    2   10
[3,]    3   20
[4,]    4   10
[5,]    5   10
[6,]    6   20

Then reclassify:

> as.matrix(reclassify(x,map))
     [,1] [,2] [,3] [,4] [,5]
[1,]   10   10   20   10   10
[2,]   20   20   20   20   20
[3,]   20   20   20   20   20
[4,]   30   30   30   30   30
[5,]   30   30   30   30   30

Note I've converted to matrix to show the change - you can just assign to a new object. You can also put NA in the second column for a NODATA value.

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.