I have MODIS images for the EVII index, with a polygon shapefile for specific areas in my area :



all polygons

all polygons



The images of NDVI and the polygons have the same projection:

> proj4string(data.shape)
[1] "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
> proj4string(r)
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

My problem is that I can't overlay my polygon on my NDVI images.

I tried this script but did not work:


r <- raster("D:/Rteledetection/Pivots/MODIS/MOD13Q1.006__250m_16_days_EVI_doy2000209_aid0001.tif")
str_name <- "MOD13Q1.006__250m_16_days_EVI_doy2000209_aid0001.tif"
imported_raster <- raster(str_name)

# importer les fichiers shapefil
data.shape <- readOGR(dsn = "D:/Rteledetection/Pivots/", layer = "2001")
plot(data.shape, border = "red")

# vector-raster
plot(data.shape, col = "red", border = "black", add = TRUE)

I would like to know how I can visualize my polygon above the image?

I have already visualized my image and the polygons above on the image in QGIS and it works but I have to zoom on the polygons to visualize them well, the thing that does not exist can not be on R.

enter image description here

enter image description here enter image description here

@BrunoConteLeite I tried your suggestion and it works, but i have another question. Do I have to load the new line for projection every time? - another question is that when I tried your script to zoom on the limit of the polygon it gives me results like :

enter image description here

that, even if for the 1st polygon when I zoom it goes beyond the limits:

enter image description here

For my objective I must extract the spectral profiles inside my polygon , but I do not know if I will have to delimit the area where are the polygons or just zoom before extract the spectral profiles?

  • 1
    It will be easier to help if you can provide sample of your data
    – Tung
    Nov 26 '18 at 7:11
  • If you cant provide sample data then at least show the output of "summary" on the data, or anything that gives us the bounding box and projections of the objects. Otherwise we can't really help.
    – Spacedman
    Nov 26 '18 at 8:20
  • yes i can provide my data : i have 2 data 1 for imag MODIS EVI and second for shapfil polygon for some areas in my area
    – Tiera Eyek
    Nov 26 '18 at 11:20
  • my problem i can't visualize my plygon above the imag , may be because it so small ? sorry i can't speak english verry well
    – Tiera Eyek
    Nov 26 '18 at 11:32
  • @Spacedman how can i do this
    – Tiera Eyek
    Nov 26 '18 at 12:00

Suppose you want to plot the raster.tif file and the polygons contained in the MLI_adm0.shp shapefile of this hypothetical data. Load them in R


r <- raster('raster.tif')
data.shape <- readOGR('.', 'MLI_adm0')

As you mentioned, you must make sure that they have the same projection. It seems that it is not the case in your data -- your data.shape seems to be projected at UTM zone 21 and your raster not (double check that). To make sure they have the same projection, use projectRaster() for the raster and spTransform() for the polygon(s) as follows

r <- projectRaster(r, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84")
data.shape <- spTransform(data.shape, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

make sure that they are correctly projected now

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

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

Now they can be put in a plot together. Use the options xlim and ylim to "zoom into" the area of your polygons. If not, the area in the plot is the one of the raster, as in the figure below.

plot(data.shape, add = T)

plot with no zoom

A nice trick to zoom the plot within the boundaries of the polygon data.shape is to use the bbox() function to obtain the boundaries automatically. Set the limits when calling the plot of the raster, as follows

    xlim=c(bbox(data.shape)[1], bbox(data.shape)[3]),
    ylim= c(bbox(data.shape)[2], bbox(data.shape)[4])
plot(data.shape, add = T)

You should obtain something like

plot with zoom

  • i have posted i result , but i have some question
    – Tiera Eyek
    Nov 27 '18 at 0:19

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.