[GIS] R: shapefile extent larger than raster extent

coordinate systemextentsrraster

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 compareCRS() gives 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[1], 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)

Best Answer

No, the shapefile extent isn't bigger than the whole world.

The problem here is the raster file, let's check the summary:

> 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

Extent from 0 to 7200 and 0 to 3600 in a UTM projection? This means that your raster is 7.2 Km x 3.6 Km, too small for your description.

Is a MODIS file (indeed 0.05 x 0.05 degrees)? MOD13 C1, MOD13 C2, MYD13 C1 or MYD13 C2? I would change extent to -180, 180, -90, 90 and CRS to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs. Try this and reproject your shapefile to the same CRS.

Why do I say this? Check product description (from MODIS products mentioned above):

Temporal Coverage                   February 18, 2000 -
Area                                Global
File Size                           ~35.5 MB
Projection                          Lat/Lon
Data Format                         HDF-EOS
Dimensions                          3600 x 7200 rows/columns
Resolution                          0.05 degree (5600-meter)
Science Data Sets (SDS HDF Layers)  13
Related Question