problem with extract() raster package using weight=TRUE

classic Classic list List threaded Threaded
8 messages Options
up
Reply | Threaded
Open this post in threaded view
|

problem with extract() raster package using weight=TRUE

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:

> t500 <-raster("t500.nc")
> projection(t500)
[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

> regioni <- readOGR(".", "ITA_adm1")
#ITA_adm1 from GADM database of Global Administrative Areas, shapefile format

> summary(regioni)
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]
(...)

> lombardia <-regioni[regioni$NAME_1=="Lombardia",]
> liguria <-regioni[regioni$NAME_1=="Liguria",]

Now, I extract values from netcdf file that fall into lombardia:

>extract(t500,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:

> extract(t500,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.

> extract(t500,lombardia,weight=T)
[[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:

> extract(t500,liguria,weight=T)
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
Reply | Threaded
Open this post in threaded view
|

Re: problem with extract() raster package using weight=TRUE

Robert 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 <[hidden email]> wrote:

> 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:
>
>> t500 <-raster("t500.nc")
>> projection(t500)
> [1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
>
>> regioni <- readOGR(".", "ITA_adm1")
> #ITA_adm1 from GADM database of Global Administrative Areas, shapefile
> format
>
>> summary(regioni)
> 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]
> (...)
>
>> lombardia <-regioni[regioni$NAME_1=="Lombardia",]
>> liguria <-regioni[regioni$NAME_1=="Liguria",]
>
> Now, I extract values from netcdf file that fall into lombardia:
>
>>extract(t500,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:
>
>> extract(t500,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.
>
>> extract(t500,lombardia,weight=T)
> [[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:
>
>> extract(t500,liguria,weight=T)
> 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.
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
up
Reply | Threaded
Open this post in threaded view
|

Re: problem with extract() raster package using weight=TRUE

up
Robert Hijmans wrote
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
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
Reply | Threaded
Open this post in threaded view
|

Re: problem with extract() raster package using weight=TRUE

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
Assistant Profesor in Applied Economics and Management
Cornell University
Reply | Threaded
Open this post in threaded view
|

Re: problem with extract() raster package using weight=TRUE

Robert 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
<[hidden email]> wrote:

> 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.
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
up
Reply | Threaded
Open this post in threaded view
|

Re: problem with extract() raster package using weight=TRUE

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
Reply | Threaded
Open this post in threaded view
|

Re: problem with extract() raster package using weight=TRUE

Ariel Ortiz-Bobea
This post was updated on .
In reply to this post by Robert Hijmans
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 helps 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
Assistant Profesor in Applied Economics and Management
Cornell University
up
Reply | Threaded
Open this post in threaded view
|

Re: problem with extract() raster package using weight=TRUE

up
Ariel Ortiz-Bobea wrote
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 helps others looking to do the same. See code below:
(...)
Many thanks Ariel for your suggestion! :-)
I'm going to use your piece of code to my purposes.
Buona giornata!
Umberto