Shapefile and Raster with same CRS but NO output when clipping

Question

Problem: Plot is empty, no results when clipping raster with shapefile.

I am using a shapefile with original EPSG4326 projection and a MODIS product with original sinusoidal projection. I converted both to the same projection (DesiredCRS) as you can see in the script, however when making a clip of the raster I don't get any results.

library(terra)
library(sf) 
library(sp) 
library(raster)

# Inputs
HDFfile <- "ModisProductsOriginal/MCD18A1.A2001001.h15v05.061.2020097222704.hdf"
Shapefile <- "Shapefile/Outline_5/Portugal_Outline_5_CAOP2019.shp"
DesiredCRS <- "+proj=longlat +datum=WGS84"

# Feature
Shp <- st_read(Shapefile)
Terceira <- Shp[Shp$DI == '43',1]
Terceira <- st_transform(Terceira, DesiredCRS)
plot(Terceira, axes=TRUE)

Terceira's plot: Plot Output

# Modis data
SDSs <- sds(HDFfile)
SDS8 <- SDSs[8]
SDS8_template <- rast(ncol=1200, nrow=1200, xmin=-34, xmax=-24, ymin=36, ymax=41, crs=DesiredCRS)
SDS8_reprojected <- project(SDS8, SDS8_template) # Reproject changes pixels?
SDS8_raster <- raster(SDS8_reprojected)
plot(SDS8_raster)

SDS8_raster's plot: Plot Output

# Clip
Clip.step1 <- crop(SDS8_raster, extent(Terceira))
Clip <- mask(Clip.step1, Terceira)
plot(Clip)

Clip.step1's plot: Plot Output

Clip's plot: Plot Output

All images show axes with the same projection. I don't understand what I'm doing wrong... Am I correctly defining the CRS and the extent?

Updated:

Panoply view of original shapefile and original HDF product: Panoply Output

Plot of the two datasets together by Hijmans: Plot Output

Answer 1

This is but tricky to debug because you mix different packages, and you do not show(object). Anyway, here is a terra approach, showing where it would be useful to see the objects metadata; and a plot that shows the raster and vector data together.

library(terra)
HDFfile <- "MCD18A1.A2001001.h15v05.061.2020097222704.hdf"
Shapefile <- "Portugal_Outline_5_CAOP2019.shp"
DesiredCRS <- "+proj=longlat +datum=WGS84"

Shp <- vect(Shapefile)
Terceira <- Shp[Shp$DI == '43',1]

SDSs <- sds(HDFfile)
SDS8 <- SDSs[8]

Start by projecting the vector data to the original raster data.

TercSin <- project(Terceira, crs(SDS8))
plot(SDS8, 1)
lines(TercSin)

As you saw, that does not work. The reason is that GDAL/PROJ reads from file or assumes the wrong ellipsoid

cat(substr(crs(SDS8),1,230), "\n")
#PROJCRS["unnamed",
#    BASEGEOGCRS["Unknown datum based upon the Clarke 1866 ellipsoid",
#        DATUM["Not specified (based on Clarke 1866 spheroid)",
#            ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
               

Or like this

describe("HDF4_EOS:EOS_GRID:\"MCD18A1.A2001001.h15v05.061.2020097222704.hdf\":MODISRAD:GMT_1200_DSR")[1:8]

#[1] "Driver: HDF4Image/HDF4 Dataset"                                         
#[2] "Files: MCD18A1.A2001001.h15v05.061.2020097222704.hdf"                   
#[3] "Size is 1200, 1200"                                                     
#[4] "Coordinate System is:"                                                  
#[5] "PROJCRS[\"unnamed\","                                                   
#[6] "    BASEGEOGCRS[\"Unknown datum based upon the Clarke 1866 ellipsoid\","
#[7] "        DATUM[\"Not specified (based on Clarke 1866 spheroid)\","       
#[8] "            ELLIPSOID[\"Clarke 1866\",6378206.4,294.978698213898,"      

As you point out (and see here), MODIS (or at least some products) uses

modcrs <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m"

So let's try that

TercSin <- project(Terceira, modcrs)
plot(SDS8, 1)
lines(TercSin)

That looks OK. So we need to do

crs(SDS8) <- modcrs

and continue

Terceira <- project(Terceira, DesiredCRS)
Terceira

SDS8_template <- rast(ncol=1200, nrow=1200, xmin=-34, xmax=-24, ymin=36, ymax=41, crs=DesiredCRS)
SDS8_reprojected <- project(SDS8, SDS8_template) 
SDS8_reprojected

# Again plot the two datasets together
plot(SDS8_reprojected)
lines(Terceira)

Clip.step1 <- crop(SDS8_reprojected, Terceira)
Clip <- mask(Clip.step1, Terceira)
Clip

plot(Clip)
lines(Terceira)