0

I have hundreds of shapefiles without a coordinate reference system. My goal is the overlay the spatial polygons over the WorldClim raster layer. I used this approach before without any problems. However, this time the coordinates from my shapefiles are strange for me. Each coordinate for bounding box and coords within polygons is composed of 8 digit numbers without comma or dot to separate de decimals. This is the bounding box from one of the shapes: 

SHP bbox: xmin:-17367529, xmax:17367529, ymin:-5997367 and ymax:7052489

which are clearly different from the bounding box of the WorldClim raster layer.

WorldClim bbox: xmin=-180,xmax=180,ymin=-60 and ymax=90

When I tried to overlay the shapefile over the raster layer using plot command nothing happens.

plot(shapefile, add=T)

I understood that this is a projection problem. Then I tried to assign the same coordinate system of the WorldClim raster layer in the shapefile using the CRS function. However, the result remains the same (i.e. the shapefiles do not over the raster). In the sequence, I tried to use the spTransform function from the rgdal package to reproject the shapefile coordinates. However, because shapefile does not have any reference system the function does not work and I do not know how to reproject the shapefile in order to match with the raster layer. I've been researching for a few days about how to deal with this problem and I believe that the absence of a reference system is a key point to the problem. However, I'm failing to overcome this problem and I would like to know if someone could help how to deal with this situation.

Ricardo Adelino
  • 165
  • 1
  • 1
  • 12

1 Answers1

2

You need to define the projection of shape files first using proj4string(meuse) or crs(shapefile)<-crs string then you can use spTransform:

library(rgdal)
data(meuse)
coordinates(meuse) <- c("x", "y")

Here you have the spatial data with x and y but you do not have the crs yet! So if you use spTransform it will fail. 

summary(meuse) #proj4string : [NA] so below line fails!
meuse.utm <- spTransform(meuse, CRS("+proj=utm +zone=32 +datum=WGS84"))

# Error in spTransform(xSP, CRSobj, ...) : 
#   No transformation possible from NA reference system

To get around this, as mentioned above, you first need to define the projection as below:

proj4string(meuse) <- CRS(paste("+init=epsg:28992",
"+towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812"))
summary(meuse) #proj4string : epsg:28992...  and then you may use spTransform

and then:

meuse.utm <- spTransform(meuse, CRS("+proj=utm +zone=32 +datum=WGS84"))
Majid
  • 1,836
  • 9
  • 19
  • spTranform approach does not work on defining projection. I believe that I'm failing because the coordinates of shapefiles are not in lonlat. The WorldClim website indicates that the available rasters do not have projections and latitude and longitude _lonlat_ coordinates were used to define the raster CRS. According to the spatial reference website [link](https://spatialreference.org/ref/epsg/wgs-84/) the longlat is equivalent to EPGS4326. I've tried to define the projection of shapefile for EPGS4326 without success. – Ricardo Adelino Nov 11 '19 at 18:42
  • I did not say: *spTranform approach may be used to define projection.* what I said is you need to define the original `crs` for your shapefile (which can be something else rather than lonlat). `EPGS4326` works for WorldCom but if it doesn't work for your shapefile then it does not relate to it! You have to find the original `crs` for your shapefile which is not really complicated. Define it to **what it is** then you can **transform** it!! Defining and transforming are not the same (look at my *example* script). – Majid Nov 11 '19 at 23:20