Sorry for top-posting, travelling: Set the environment variable PROJ_NETWORK=ON, and start R, or Sys.setenv("PROJ_NETWORK"="ON") before loading sf:
library(sf) station_pts<-read.csv(text="Station_ID, Lat, Lon 1, 55.8533, -2.38588 2, 55.8551, -2.39620 3, 55.8473, -2.11304") stations<-st_as_sf(station_pts, coords=c("Lon", "Lat"), crs="OGC:CRS84") sf_proj_network() sf_proj_network(FALSE) o <- sf_proj_pipelines("OGC:CRS84", "EPSG:27700") o$definition o$accuracy stations_OS_no_cdn <- st_transform(stations, "EPSG:27700", pipeline=o$definition[1]) stations_OS_no_cdn sf_proj_network(TRUE) stations_OS <- st_transform(stations, "EPSG:27700") stations_OS There is about a 1m difference between the transformation using a grid and the best Helmert transformation. My guess is that even after other corrections, the online point-at-a-time transformation pages use a grid rather than a Helmert transformation. There is usually a reason for observed differences, and https://cdn.proj.org is a very useful resource where mapping agencies have contributed grids. Roger -- Roger Bivand Emeritus Professor Norwegian School of Economics Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway roger.biv...@nhh.no ________________________________________ From: R-sig-Geo <r-sig-geo-boun...@r-project.org> on behalf of Micha Silver <tsvi...@gmail.com> Sent: 23 October 2022 13:58 To: Nick Wray; r-sig-geo@r-project.org Subject: Re: [R-sig-Geo] Converting longitude and latitude to UK OS references Just to combine Rob's and Roger's responses: library(sf) # Create a data.frame from text station_pts<-read.csv(text="Station_ID, Lat, Lon + 1, 55.8533, -2.38588 + 2, 55.8551, -2.39620 + 3, 55.8473, -2.11304") # Make an sf object from the DF # Careful about the order of Lon/Lat stations<-st_as_sf(station_pts, coords=c("Lon", "Lat"), crs=4326) # CRS transform stations_OS <- st_transform(stations, 27700) str(stations_OS) Classes ‘sf’ and 'data.frame': 3 obs. of 2 variables: $ Station_ID: int 1 2 3 $ geometry :sfc_POINT of length 3; first list element: 'XY' num 375940 662301 - attr(*, "sf_column")= chr "geometry" - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA ..- attr(*, "names")= chr "Station_ID" # Just to be sure, display CRS: st_crs(stations_OS) Coordinate Reference System: User input: EPSG:27700 wkt: PROJCRS["OSGB 1936 / British National Grid", BASEGEOGCRS["OSGB 1936", DATUM["OSGB 1936", ELLIPSOID["Airy 1830",6377563.396,299.3249646, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4277]], CONVERSION["British National Grid", METHOD["Transverse Mercator", ID["EPSG",9807]], PARAMETER["Latitude of natural origin",49, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8801]], PARAMETER["Longitude of natural origin",-2, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8802]], PARAMETER["Scale factor at natural origin",0.9996012717, SCALEUNIT["unity",1], ID["EPSG",8805]], PARAMETER["False easting",400000, LENGTHUNIT["metre",1], ID["EPSG",8806]], PARAMETER["False northing",-100000, LENGTHUNIT["metre",1], ID["EPSG",8807]]], CS[Cartesian,2], AXIS["(E)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["(N)",north, ORDER[2], LENGTHUNIT["metre",1]], USAGE[ SCOPE["Engineering survey, topographic mapping."], AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."], BBOX[49.75,-9,61.01,2.01]], ID["EPSG",27700]] If you have your station Long/Lat coordinates in a CSV file, then `st_read()` will read directly into an sf object. On 23/10/2022 12:08, Nick Wray wrote: > I am trying to convert longitude and latitude values for UK weather > stations to UK Ordnance Survey National Grid References. There are sites > where one can do them one at a time but I have a large number. I have > found some code which does the conversion and include the first three > points which I want to convert as an example > > library(sp) > xy<-as.data.frame(cbind(c(55.8533,55.8551,55.8473),c(-2.38588,-2.39620,-2.11304))) > colnames(xy)<-c("lon","lat") > xy > coordinates(xy)<-~lon+lat > ## see site in jounral convert llon lat > proj4string(xy)<-CRS("+init=epsg:4326") > ptsOS<-spTransform(xy,CRS("+init=epsg:27700")) > ptsOS > > but it doesn't give the right answers - the output doesn't correspond to > either what I get whn I put an individual point into a conversion site nor > with other points in the same region which I have got from other sources > and are mutually compatible > > Can anyone explain where I'm going wrong and how I can fix this? Thanks > Nick Wray > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&data=05%7C01%7CRoger.Bivand%40nhh.no%7C68be57c5516345d3e61a08dab4edf0ee%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638021231300728219%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=mqCLcKYV2K%2FuxI3wmkhC6QMbALy5Y7uuVXBHkVoyAxE%3D&reserved=0 -- Micha Silver Ben Gurion Univ. Sde Boker, Remote Sensing Lab cell: +972-523-665918 https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Forcid.org%2F0000-0002-1128-1325&data=05%7C01%7CRoger.Bivand%40nhh.no%7C68be57c5516345d3e61a08dab4edf0ee%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638021231300728219%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=5IJ6yY1pJaYtRglQLNV%2FJY3SfktRn0864FPres66CSI%3D&reserved=0 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&data=05%7C01%7CRoger.Bivand%40nhh.no%7C68be57c5516345d3e61a08dab4edf0ee%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638021231300728219%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=mqCLcKYV2K%2FuxI3wmkhC6QMbALy5Y7uuVXBHkVoyAxE%3D&reserved=0 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo