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

lutz Published at Dev

Lutz

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?

Chris

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')

Collected from the Internet

Please contact [email protected] to delete if infringement.

edited at2021-05-26

Comments

0 comments

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

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

pump.io port in URL

How to import an asset in swift using Bundle.main.path() in a react-native native module

Failed to listen on localhost:8000 (reason: Cannot assign requested address)

Double spacing in rmarkdown pdf

SQL Server : need add a dot before two last character

C++ 16 bit grayscale gradient image from 2D array

JMeter: Why get error when try to save test plan

JWT gives JsonWebTokenError "invalid token"

How to make thrown errors visible outside of a Promise?

How to tell if iOS Today Widget is being updated in the background?

Calling Doctrine clear() with an argument is deprecated

Capybara Selenium Chrome opens About Google Chrome

How to update azerothcore-wotlk docker container

Adding Ripple Effect to RecyclerView item

mysql.connector.errors.InterfaceError: 2003: Can't connect to MySQL server on '127.0.0.1:3306' (111 Connection refused)

Error while applying filter on dataframe - PySpark

Unable to add slack to bluemix project

MyPy fails dataclass argument with optional list of objects type

How can I validate and parse phone numbers to extract their country calling code and area code?

Single Sign-On in Spring by using SAML Extension and Shibboleth

python how to create many-to-many of lists inside one list