I am using R for my spatial analysis. I am fairly new to this so I am still trying to find my way around things like changing projections etc...
I have a raster with vegetation index values (0.5 x 0.5 degrees), and a shapefile with polygons. I wanted to use
intersect() or something along those lines to get the vegetation index values for my polygons, but I am running into issues with the extent or projection - not exactly sure what is causing the issue.
At first it became clear the .tif didn't have a projection, so I set the projection the same as my shapefile. Testing with
TRUE. When I try to do
intersect() though, I get the error
Warning messages: 1: In intersect(x, y) : non identical CRS 2: In intersect(x, y) : polygons do not intersect
When inspecting the raster and shapefiles, I see indeed that the extents are different - but the extent of my shapefile is actually larger than my raster (and my raster covers the whole world, I plotted it and that's correct). The shapefile is a region in the Amazon, so it should fall inside the raster covering the whole world.
How does that happen? Can I change the extent of my raster without messing up the resolution or the grids?
I've done a lot of searching online, but most of the resources talk about cropping when the extent of one file falls within the other file...
Here are the details of my files (
r is the raster,
comms the shapefile):
> show(r) class : RasterLayer dimensions : 3600, 7200, 25920000 (nrow, ncol, ncell) resolution : 1, 1 (x, y) extent : 0, 7200, 0, 3600 (xmin, xmax, ymin, ymax) coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 data source : /EVI2_testdata/EVI2_1983152.tif names : EVI2_1983152 values : -32768, 32767 (min, max) > show(comms) class : SpatialPolygonsDataFrame features : 99 extent : 320877.6, 703893, 8549725, 8916662 (xmin, xmax, ymin, ymax) coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 variables : 9 names : Join_Count, TARGET_FID, P5_1, TabulateAr, Shape_Leng, Shape_Area, Cluster_ID, ClusterID1, Cluster1v2 min values : 0, 0, 3 Arroyos, 1, 276.5641, 1.010235e+08, 1, 1, 0 max values : 19, 4282, Virgen de la Candelaria, 103, 217329.5623, 9.944712e+07, 4, 4, 0
Below the code used to get to this point:
library(gdalUtils) library(raster) library(rgdal) # Get a list of sds names out of HDF4 file, make a tif sds <- get_subdatasets('VIP30P2.A1983152.004.2014270210035.hdf') filename <- 'EVI2_1983152.tif' gdal_translate(sds, dst_dataset = filename) # Load the Geotiff created into R r <- raster(filename) # read in shapefiles comms<-readOGR(dsn=path.expand("~/MAP_PhD_research/DATA/Shapes/CommunityPolygons_new/MAP_Tenure_WithComminutyPnt_20130709.shp"), layer="MAP_Tenure_WithComminutyPnt_20130709") # check projections identicalCRS(comms,r) # which is false proj4string(comms) # result: +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 # assign same projection to raster file crs(r)<-"+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" proj4string(r) compareCRS(r,comms) # result is true # still intersect does not work intersect(comms,r)