1

It's a shapefile that has a column with population.

a$POP2001

It has 200 rows.I want to extract the settlements that contain less than 1000 residents, so i did this:

s<-a[a$POP2001<1000,]

Even though it has 200 rows, when i run a$POPULATION it returns thousands of rows. After looking carefully i saw that it duplicates multiple times the 200 rows.So by plotting the "s" it shows all the settlements, render it useless to extract numbers. What should i do to make it return the actual settlements that are in accordance with the previous statement?

str(a)


Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
 ..@ data       :'data.frame':  25413 obs. of  24 variables:
 .. ..$ ZX        : Factor w/ 2 levels "83160101","83160701": NA NA NA NA NA        NA NA NA NA NA ...
 .. ..$ OBJECTID  : num [1:25413] 1116 1116 1116 1116 1116 ...
 .. ..$ CODE_OIK  : Factor w/ 199 levels "83010101","83010102",..: 182 182    182 182 182 182 182 182 182 182 ...
 .. ..$ NAME_OIK  : Factor w/ 193 levels "ΠαλαιοχΟΟΞΉΞΏΞ½""|   __truncated__,..: 111 111 111 111 111 111 111 111 111 111 ...
 .. ..$ CODE_GDIAM: Factor w/ 1 level "8": 1 1 1 1 1 1 1 1 1 1 ...
 .. ..$ NAMEF_OIK : Factor w/ 194 levels "ΠαλαιοχΟΟΞΉΞΏΞ½,το""|   __truncated__,..: 111 111 111 111 111 111 111 111 111 111 ...
  .. ..$ POINT_X   : num [1:25413] 708229 708229 708229 708229 708229 ...
  .. ..$ POINT_Y   : num [1:25413] 4317433 4317433 4317433 4317433 4317433   ...
  .. ..$ LAT       : num [1:25413] 39 39 39 39 39 ...
  .. ..$ LON       : num [1:25413] 26.4 26.4 26.4 26.4 26.4 ...
  .. ..$ H         : num [1:25413] 80 80 80 80 80 80 80 80 80 80 ...
  .. ..$ EDRA_DIAM : num [1:25413] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ CODE_DIAM : Factor w/ 73 levels "83010100","83010200",..: 68 68 68   68 68 68 68 68 68 68 ...
  .. ..$ NAME_DIAM : Factor w/ 73 levels "Τ.Δ.Παλαιοκήπου",..:   49 49 49 49 49 49 49 49 49 49 ...
  .. ..$ CODE_OTA  : Factor w/ 13 levels "83010000","83020000",..: 12 12 12   12 12 12 12 12 12 12 ...
  .. ..$ NAME_OTA  : Factor w/ 13 levels "ΔΗΞΞΞ£ ΠΛΩΞΑΡΙΞΞ¥""|   __truncated__,..: 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ CODE_NOM  : Factor w/ 1 level "83": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ NAME_NOM  : Factor w/ 1 level "ΞΞΞΞΞ£ ΛΕΣΒΞΞ¥""| __truncated__: 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ NAME_GDIAM: Factor w/ 1 level "ΞΞ—Ξ£ΞΞ™ ΑΙΓΑΙΞΞ¥""| __truncated__: 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ DIST1     : int [1:25413] 500 500 500 500 500 500 500 500 500 500 ...
  .. ..$ DIST2     : int [1:25413] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 ...
  .. ..$ CODE2001  : num [1:25413] 8.3e+07 8.3e+07 8.3e+07 8.3e+07 8.3e+07 ...
  .. ..$ POP2001   : num [1:25413] 192 172 426 177 97 76 4 53 465 465 ...
  .. ..$ vat       : num [1:25413] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:25413, 1:2] 708229 708229 708229 708229 708229 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:25413] "2" "2" "2" "2" ...
  .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  ..@ bbox       : num [1:2, 1:2] 658196 4315330 726037 4360718
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr NA

length(unique(a$POP2001))
[1] 152

length(unique(paste(a$POINT_x,a$POINT_Y)))
[1] 198

length(unique(a$POP2001[which(paste(a$POINT_x,a$POINT_Y)=="7082294317433")]))
[1] 0


library(foreign)
oikismoi<-readShapeSpatial("Oikismoi/lsvoikismoiall")
df1 <- read.dbf("a2001pop.dbf")
m<-merge(oikismoi,df1,by.x="ZX",by.y="HJ",)
reclass = function(x){
ifelse(x < 500, 0,
     ifelse(x < 1000 & x >= 500, 5, 
            ifelse(x > 1000, 10, NA)))
}
m$vat<-as.numeric(lapply(m$POP2001,reclass))
qw<-m[m$POP2001>5000,]

The problem with the duplicates happened after the merge. Before it worked fine.When the shp got the POP2001.

NEW**********************************************

How i created the shapefile. Its a shapefile and a dbf. The dbf has the pop2001.I did a merge of the attributes of the shp and the dbf through a common column so the shp can contain the POP2001. Before the merge the <1000 worked well in the initial data.

   df1 <- read.dbf("a2001pop.dbf") 
   a<-merge(oikismoi,df1,by.x="ZX",by.y="HJ")
   df1[df1$POP2001>5000,]
   CODE2001 POP2001       HJ
   1   83010000   37881 83160101
   2   83010100   30576 83160701
   3   83010101   28950     <NA>
   67  83050000    6945     <NA>
   96  83060000    5111     <NA>
   150 83080000    8073     <NA>

After the merge this

a$POP2001
 [1] 37881 30576 28950  6945  5111  8073 28950  6945  5111  8073 28950  6945   5111  8073
 [15] 28950  6945  5111  8073 28950  6945  5111  8073 28950  6945  5111  8073 28950  6945
 [29]  5111  8073 28950  6945  5111  8073 28950  6945  5111  8073 28950  6945  5111  8073
 [43] 28950  6945  5111  8073 28950  6945  5111  8073 28950  6945  5111  8073 28950  6945 AND SO THE LIST GOES ON...
  • can you post what you get from str(a)? – aaryno Jul 13 '15 at 23:31
  • 2
    You HAVE got to provide reproducible examples or in this case enough previous code for us to see where things may have gone wrong. – Jeffrey Evans Jul 14 '15 at 1:43
  • i tried it on the meuse and it worked well.Any ideas what's going on? – gsa Jul 14 '15 at 9:24
  • I'm assuming you mean a$POP2001 because there's no a$POPULATION. I always want to point out that this data source you are using already has the duplicated data. I'm assuming that the fields POINT_X and POINT_Y represent coordinates. It is apparent that the first 5 values are identical, yet the POP2001 values are different. I'm going to guess that the data has become scrambled before you got to this point. Tell us what it says for the following commands: – aaryno Jul 14 '15 at 16:47
  • length(unique(a$POP2001)) and length(unique(paste(a$POINT_x,a$POINT_Y))) and length(unique(a$POP2001[which(paste(a$POINT_x,a$POINT_Y)=="708229 4317433")])). If my suspicions are correct you will have about 199 unique values in all three cases. How did the data get created? – aaryno Jul 14 '15 at 16:47
2

Merge is not doing what you want, probably because there are only two unique values in HJ or ZX. The merge columns should uniquely identify each record in oikismo and df1. Change your by.x and by.y appropriately.

  • i just added manually the POP2001 to the shapefile's attibute table and it worked great. Thanks – gsa Jul 15 '15 at 10:08
1

You should note that str is telling us that you have 25413 observations, not 200, in your "a" @data (data.frame) object. You could have a multi-part shapefile where numerious records are associated with each spatial location. If the results of dim(a) and dim(a@data) differ this is likely the case. This is why you provide as much information as possible in you initial post.

You also have corrupt or special character data in "NAME_OIK", "NAMEF_OIK", "NAME_DIAM", "NAME_OTA", "NAME_NOM" and "NAME_GDIAM". This could cause you issues in the future. You may also want to use "readOGR" in "rgdal" rather then the maptools package.

Let's step through an analysis that mirrors what you are attempting. Please note that readShapeSpatial will result in an sp object and read.dbf will result in a data.frame object so, the objects in my example will represent the same object classes. Read over syntax very carefully and take note of the S4 slot format with the usage of "@" to index objects in various slots of the sp object.

Create a SpatialPointsDataFrame object "a" with a POP2001 column (n=1-10000)

library(sp)
data(meuse)
a <- meuse
coordinates(a) <- ~x+y
a@data <- data.frame(ZX=1:nrow(a), POP2001=round(runif(nrow(a), 1, 10000),0) )

Classify the POP2001 column and add results "vat" to @data slot in "a"

a@data <- data.frame(a@data, vat = ifelse(a@data$POP2001 < 500, 0, 
          ifelse(a@data$POP2001 >= 500 & a@data$POP2001 < 1000, 5, 
          ifelse(a@data$POP2001 >= 1000, 10, NA))))

Create a new data.frame "df1", merge with "a" and display first 6 rows (I don't know what this has to do with your problem but, here it is).

df1 <- data.frame(HJ=a@data$ZX, y=runif(nrow(a)) )
a <- merge(a, df1, by.x="ZX", by.y="HJ")
head(a@data)

Subset "a" to a new "lt1000" object representing POP2001 < 1000

lt1000 <- a[a@data$POP2001 < 1000 ,]  
dim(lt1000)
dim(a)

Plot results

plot(a, pch=20, col="black")
  plot(lt1000, pch=20, col="red", add=TRUE)
  box()
  legend("topleft", legend=c(">= 1000", "< 1000"),
        pch=c(20,20), col=c("black","red")) 
  • as i said the merge part happerns so the shp can have the POP2001 from the dbf.This is the part where the thing goes bad.Do u know whats going on? – gsa Jul 14 '15 at 23:10
  • First, you have an extra comma in your merge syntax. Second, make sure that you have the sp library added. As a test, I have added a bunch of random rows to "df1" and still got a one-to-one match. Can you run through the worked example without issues? It would help if you posted your data. – Jeffrey Evans Jul 14 '15 at 23:13
  • Definitely. how can i post the data?Through email? – gsa Jul 14 '15 at 23:17
  • What you did in the answer is the same as the question but with different data – gsa Jul 14 '15 at 23:28

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.