1

I have two set of data:

  1. raster land cover from [Globcover]:(http://www.fao.org/land-water/land/land-governance/land-resources-planning-toolbox/category/details/en/c/1036356/)
  2. A set of point shapefiles of health facilities in Africa.[Google drive]:(https://drive.google.com/drive/folders/1Oi9Bd_T_W82DXZjN7XUCoAyvM6kCVkYQ?usp=sharing)

I have problems with ensuring that they are aligned after transforming their CRS from geographic coordinate system to projected coordinate system.

Both files has CRS of WGS84. I tried to plot them both in WGS84 and they are aligned perfectly, as shown below.

crs(landcover)

CRS arguments:
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

crs(healthfac_point)

CRS arguments:
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

enter image description here

However, when I convert both the CRS of the globcover data and point shapefiles into epsg:32735, they no longer line up.

healthfac_point_projected <- spTransform(healthfac_point, CRS("+init=epsg:32735"))

crs(healthfac_point_projected)

CRS arguments:
+init=epsg:32735 +proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0 

projection(landcover) <- CRS("+init=epsg:32735")

crs(landcover)
CRS arguments:
+init=epsg:32735 +proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0 

extent(healthfac_point_projected)

class       : Extent 
xmin        : -5832574 
xmax        : 4548844 
ymin        : 6136173 
ymax        : 13175174 

extent(landcover)

class       : Extent 
xmin        : -180.0014 
xmax        : 179.9986 
ymin        : -64.99861 
ymax        : 90.00139  

I also tried this.

healthfac_point_projected <- spTransform(healthfac_point, CRS("+init=epsg:32735"))

crs(healthfac_point_projected)

CRS arguments:
+init=epsg:32735 +proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

landcover_proj <- projectRaster(landcover, proj4string(healthfac_point_projected))

Error in (function (classes, fdef, mtable)  : 
unable to find an inherited method for function ‘res’ for signature ‘"character"’
In addition: Warning message:
In min(dim(to)[1:2]) : no non-missing arguments to min; returning Inf

Ultimately I would like to create buffers around these point shapefiles and find out the % land cover in each buffer. However, as their extent are different despite transforming their CRS, I cannot do that.

Can you point out what's wrong with the way I have transformed the CRS of both files?

2

It's this bit, is wrong it assigns the crs doesn't apply it by transforming the data:

projection(landcover) <- CRS("+init=epsg:32735")

you're expected to transform the raster i.e.

new <- raster(healthfac_point_projected)
res(new) <- c(150000, 150000) ## change resolution to suit, this is metres (the width,height of each cell)
landcoverutm <- projectRaster(landcover, new)

Untested. UTM is not a great choice and officially bad for such a wide set of longitudes, for this kind of region I'd probably use LambertAzimuthalEqualArea (or a Conic, ..., but this is not that question). Try

"+proj=laea +lon_0=40 +lat_0=-10 +datum=WGS84"

With newer raster/rgdal/sp/etc you'll want to drop the '+ellps' arg which will cause warnings and be unused.

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.