1

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)
  • r is not georeferenced, extent : 0, 7200, 0, 3600 (xmin, xmax, ymin, ymax) - this is the "index extent", i.e. a fall back when there's no information is just 0 to number of cells along each axis. Having utm on this is completely bogus, mis-applied metadata. – mdsumner Aug 6 '17 at 1:54
  • yes just look at the extents, they are clearly not aligned. You need to project the raster with the correct reference system (not necessarily the same as polygons) and then reproject it to the shapefile's CRS (or reproject shapefile to that of the raster) – Sam Aug 10 '17 at 8:54
1

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
  • Thanks, it is a MODIS derived product (adapted by University of Arizona). I'll give this a shot and update – Geraldine Sep 13 '17 at 17:05

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.