Unknown CRS in QGIS when projecting to EPSG:25833 in R

Question

I want to project a spatial data frame to EPSG 25833 in R but QGIS does not seem to know it (for reproducibility, I use the code jazzurro created in his/her answer to this question with minor changes)

library(rgdal)

mydf <- structure(list(longitude = c(128.6979, 153.0046, 104.3261, 124.9019, 
                                     126.7328, 153.2439, 142.8673, 152.689), latitude = c(-7.4197, 
                                                                                          -4.7089, -6.7541, 4.7817, 2.1643, -5.65, 23.3882, -5.571)), .Names = c("longitude", 
                                                                                                                                                                 "latitude"), class = "data.frame", row.names = c(NA, -8L))


### Get long and lat from your data.frame. Make sure that the order is in lon/lat.

xy <- mydf[,c(1,2)]


# Here I use the projection EPSG:25833
spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
                               proj4string = CRS("+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))


#Export as shapefile
writeOGR(spdf, "file location", "proj_test", driver="ESRI Shapefile",overwrite_layer = T) #now I write the subsetted network as a shapefile again

Now when I load the shapefile into QGIS it doesn´t know the projection.

Any ideas?

score 2 · Accepted Answer · answered Mar 02 '20 at 17:47

In making your SpatialPointsDataFrame:

# Wrong!
spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
                               proj4string = CRS("+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))`

You are telling the data frame what crs your points are in, so you should specify 4326 since your original data is lon/lat.

So it should be:

spdf <- SpatialPointsDataFrame(coords = xy, data = mydf,
                               proj4string = CRS("+proj=longlat +datum=WGS84"))

And then you can transform the data to another CRS using spTransform:

spTransform(spdf, CRS('+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'))

For this particular data, we get an error because one of the points doesn't convert to your target CRS:

Error in spTransform(xSP, CRSobj, ...) : failure in points 3 In addition: Warning message: In spTransform(xSP, CRSobj, ...) : 1 projected point(s) not finite

I prefer working in sf so we could also do:

library(sf)
sfdf <- st_as_sf(mydf, coords = c('longitude', 'latitude'), crs=4326, remove=F)
sfdf_25833 <- sfdf %>% st_transform(25833)

sfdf_25833
#> Simple feature collection with 8 features and 2 fields (with 1 geometry empty)
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 5589731 ymin: -19294970 xmax: 11478870 ymax: 19337710
#> epsg (SRID):    25833
#> proj4string:    +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#>   longitude latitude                   geometry
#> 1  128.6979  -7.4197 POINT (10198485 -17980649)
#> 2  153.0046  -4.7089  POINT (5636527 -19294974)
#> 3  104.3261  -6.7541                POINT EMPTY
#> 4  124.9019   4.7817  POINT (11478868 18432292)
#> 5  126.7328   2.1643  POINT (11046583 19337712)
#> 6  153.2439  -5.6500  POINT (5589731 -19158700)
#> 7  142.8673  23.3882   POINT (6353660 16093116)
#> 8  152.6890  -5.5710  POINT (5673080 -19163103)

and you can write and open with QGIS using:

write_sf(sfdf_25833, 'mysf.gpkg')

thanks for your answer! When applying your first suggestion to my data I get the same error as you. However switching to sf now would mean quite a lot of work because I started with SpatialDataFrames and some functions wouldn´t work with sf. Is there any other way to get around this error? — Lutz, Mar 03 '20 at 09:20
If you're still getting that error are you sure that 25833 is the right projection? Where abouts in the world are your real points and why are you transforming in the first place? You can also convert freely between `sf` and `sp` at any time using `st_as_sf` and `as` but note while there's no error with sf, the offending point is removed (empty) — Chris, Mar 03 '20 at 09:30
I changed to sf now and it works (its much faster anyways). Thanks for your help! — Lutz, Mar 10 '20 at 09:06

Unknown CRS in QGIS when projecting to EPSG:25833 in R

1 Answers1