Skip to content

problem with extract() raster package using weight=TRUE

8 messages · Robert J. Hijmans, Ariel Ortiz-Bobea, up

Original
up
# 
Hi everybody! I'm a new user, sorry for any kind of errors in posting :-).

I'm trying to use* extract()* from raster library to get data values from
netcdf file that are inside a shapefile.
I'm using rgdal library to manage shapefiles, and raster library to read
netcdf file (a very simple file).
This is the simple code:
[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
#ITA_adm1 from GADM database of Global Administrative Areas, shapefile
format
Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x  6.630879 18.52069
y 35.492917 47.09096
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
(...)
Now, I extract values from netcdf file that fall into lombardia:
[[1]]
  [1] 249.3597 249.3362 249.3284 249.9573 249.5706 249.5081 249.4065
250.2269
  [9] 250.2581 250.2151 250.1058 249.9964 249.9261 249.8831 249.8245
249.7347
(...)

Now into liguria:
[[1]]
 [1] 249.9729 250.3440 250.2269 250.1800 250.1331 250.0823 250.0198 249.9729
 [9] 249.9573 249.9651 250.4026 250.3128 250.2425 250.1956 250.0862 250.0550
[17] 250.0628 250.1019 250.1487 250.4808 250.4222 250.3753 250.2737 250.2972
[25] 250.3284 250.7972 250.4495 250.4456 250.4495 250.4808 250.8362 250.6800
[33] 250.5940 250.5862 251.1409 251.0003 250.8753 250.8011 250.7894

Now I try to extract with weight enabled.
[[1]]
          value weight
  [1,] 249.4651   0.23
  [2,] 249.3948   0.27
  [3,] 249.8323   0.15
  [4,] 249.8206   0.42
  [5,] 249.8206   0.10
  [6,] 249.4183   0.06
  [7,] 249.3597   0.99
  [8,] 249.3362   0.99
  [9,] 249.3284   0.90
(...)

Nice one!
Now I do the same thing on the second shapefile:
Error: identicalCRS(x, y) is not TRUE

That's the problem! I can extract from liguria netcdf values, but I cannot
extract them with weight.

Any suggestion?

Manymany thanks in advance!

Umberto



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/problem-with-extract-raster-package-using-weight-TRUE-tp7584665.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
Robert J. Hijmans
# 
Umberto,

The work-around is to assure that the CRS are the same:

projection(regioni) <- projection(t500)

You should of course only do that when the projections are in fact the
same, as in this case. The problem is caused by sp::identicalCRS that
compares the verbatim, but not the semantic representation of the CRS
(and, therefore, perhaps should give a warning rather than an error?).

Robert
On Tue, Sep 17, 2013 at 5:10 AM, up <u.p at operamail.com> wrote:
1 day later
up
# 
Robert Hijmans wrote
Many thanks Robert for your suggestion! :-)

But I've already tried that, and did not work.
Yestarday I tried to use another shapefile, because I suspected something
wrong in my WGS84 shapefile.
And with another shapefile things gone well!

I've downloaded a shapefile of my country but in a different projection
("+proj=utm +zone=32 +ellps=intl +units=m +no_defs"), and then reprojected
in "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" with:

liguria<-spTransform(ligu_utm, CRS("+proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0"))

Then I've extracted values with the same function:
estra_liguria<-extract(prec,liguria,weights=T,cellnumbers=T,df=T)

and that's the result:

   ID cell totprec weight
1   1   10       0   0.12
2   1   18       0   0.34
3   1   19       0   0.33
4   1   20       0   0.42
5   1   25       1   0.57
6   1   26       0   0.60
7   1   27       0   0.84
8   1   28       0   0.96
9   1   29       0   0.91
10  1   30       0   1.00
11  1   34      10   0.06
12  1   35       5   0.84
13  1   36       0   1.00
14  1   37       0   1.00

I don't know why the first shapefile doesn't work (I've downloaded it from
http://www.gadm.org, lonlat WGS84), but the new one does the job! :-)

In any case, thanks very much!

Umberto




--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/problem-with-extract-raster-package-using-weight-TRUE-tp7584665p7584673.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
Ariel Ortiz-Bobea
# 
Hello Umberto and Robert,

I'm actually encountering the same problem.

- when I set weights=F, the extract function works
- when I set weights=T, I get the "Error: identicalCRS(x, y) is not TRUE"
- I've used this code before with the same shapefile and raster file (with
an older version of raster) and everything worked fine. I presume the newer
versions of extract() are more sensitive
- Following the solution by Umberto, I tried several US shape files
(including the latest one from the census
http://www.census.gov/geo/maps-data/data/tiger-line.html, which should be of
good quality)... without success. I get the same error.

Any suggestions?

Thanks

Ariel

/======================
My code looks like:

# Import shape file
poly <- readOGR("tl_2013_us_county.shp", "tl_2013_us_county")
projection(poly) # I get: "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80
+towgs84=0,0,0"

# Import raster file
r   <- raster("2012010100.NOAH.grb", nl=1)
projection(r) # I get: "+proj=longlat +a=6371200 +b=6371200 +no_defs"

# change projection of polygon
poly2 <- spTransform(poly, CRS(proj4string( r )) )
projection(poly2) # As espected, I get: ""+proj=longlat +a=6371200
+b=6371200 +no_defs"

# Extract without weights (works)
 v1 <- extract( r, poly2, cellnumbers=TRUE, progress="text")

# Extract with weights (does not work)
 v2 <- extract( r, poly2, cellnumbers=TRUE, weights=T, progress="text") 
# I get: the bar advances to 1-6% depending on the shape file and then I get
"Error: identicalCRS(x, y) is not TRUE"
/







-----
Ariel Ortiz-Bobea
Fellow at Resources for the Future
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/problem-with-extract-raster-package-using-weight-TRUE-tp7584665p7584678.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
1 day later
Robert J. Hijmans
# 
The error occurs if "over" is called; whether this happens depends on
the polygons and the raster used. In earlier versions overlay was used
and there was less checking of CRS in that case. This problem has been
fixed in the next version of raster (2.1-56). For now, I think the
work around -- unlike what I said before -- is to set the CRS of both
datasets to NA.
Robert

On Thu, Sep 19, 2013 at 5:53 PM, Ariel Ortiz-Bobea
<aortizbobea at arec.umd.edu> wrote:
1 day later
up
# 
Many thanks Robert for your work, really! :-)

Another question.
With weights parameter I see the approximate fraction on of the cell is
rounded to 0.01.
Is it possible to have much resolution on that fraction, like 0.001 or more?
Perhaps is it already possible (and I've not read carefully the *man*
page...)?

Good afternoon and ciao!

Umberto



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/problem-with-extract-raster-package-using-weight-TRUE-tp7584665p7584691.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
Ariel Ortiz-Bobea
# 
That's much clever than my longer workaround using rasterize and zonal.

By the way, Umberto, my workaround could help you get more precision on the
extract weights. Hope this help others looking to do the same. See code
below:

---------------

# Create raster
   library(raster)
   r <- raster(ncol=36, nrow=18)
   r[] <- 1:ncell(r)
   cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60),
c(-180,-20))
   cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
   polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1),
         
# visualize
   plot(r)
   plot(polys, add=T)
     
# Extract the cellnumbers and weights using extract()
   v <- extract(r, polys, weights=TRUE, cellnumbers=TRUE, small=TRUE)

# Extract the cellnumbers and weights using rasterize() and zonal()
   r.mod <- r
   values(r.mod) <- 1:ncell(r)
   fac <- sqrt(1000) # this is where you set the resolution factor
   r.mod <- disaggregate(r.mod, fact=c(fac, fac)) # increases raster
resolution by a factor of "fac"
   r.polys <- rasterize(x=polys, y=r.mod) # rasterize county polygons

   # Compute weights
   w <- lapply(unique(r.polys), function(x) {
      m <- r.polys
      m[getValues(m)!=x]  <- NA
      m[getValues(m)==x]  <- 1
      out <- zonal(x=m, z=r.mod, fun="sum")
      colnames(out) <- c("cell", "weight")
      out[,"weight"] <- out[,"weight"]/(fac^2)
      out[out[,"weight"]!=0,]
     })
   names(w) <- unique(r.polys)

# See weigths
 head(v[[1]])
 head(w[[1]]) 



-----
Ariel Ortiz-Bobea
Fellow at Resources for the Future
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/problem-with-extract-raster-package-using-weight-TRUE-tp7584665p7584693.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
up
# 
Ariel Ortiz-Bobea wrote
Many thanks Ariel for your suggestion! :-)
I'm going to use your piece of code to my purposes.
Buona giornata!
Umberto



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/problem-with-extract-raster-package-using-weight-TRUE-tp7584665p7584697.html
Sent from the R-sig-geo mailing list archive at Nabble.com.