spTransform-methods.RdThe spTransform methods provide transformation between datum(s) and conversion between projections (also known as projection and/or re-projection), from one unambiguously specified coordinate reference system (CRS) to another, prior to version 1.5 using Proj4 projection arguments. From version 1.5, Well-Known Text 2 (WKT2 2019) strings are used. For simple projection, when no Proj4 +datum tags are used, datum projection does not occur. When datum transformation is required, the datum should be defined with a valid value both in the CRS of the object to be transformed, and in the target CRS. In general datum is to be prefered to ellipsoid, because the datum always fixes the ellipsoid, but the ellipsoid never fixes the datum.
In addition, before version 1.5 the +towgs84 tag should have been used where needed to make sure that datum transformation would take place. Parameters for +towgs84 were taken from the legacy bundled EPSG file if they are known unequivocally, but could be entered manually from known authorities. Not providing the appropriate +datum and +towgs84 tags led to coordinates being out by hundreds of metres. Unfortunately, there is no easy way to provide this information: the user has to know the correct metadata for the data being used, even if this can be hard to discover.
From version 1.5, spTransform uses the modern PROJ coordinate operation framework for transformations. This avoids pivoting through WGS84 if possible, and uses WKT2 (2019) strings for source and target CRS often constructed from the bundled EPSG SQLite database. The database is searched for feasible candidate coordinate operations, and the most accurate available is chosen. More details are available in a vignette: vignette("CRS_projections_transformations").
get_transform_wkt_comment()
set_transform_wkt_comment(value)
get_enforce_xy()
set_enforce_xy(value)
get_last_coordOp()default void method
returns transformed coordinates of an "SpatialPoints" object using the projection arguments in "CRSobj", of class CRS
returns transformed coordinates of an "SpatialPointsDataFrame" object using the projection arguments in "CRSobj", of class CRS
returns transformed coordinates of an "SpatialLines" object using the projection arguments in "CRSobj", of class CRS
returns transformed coordinates of an "SpatialLinesDataFrame" object using the projection arguments in "CRSobj", of class CRS
returns transformed coordinates of an "SpatialPolygons" object using the projection arguments in "CRSobj", of class CRS
returns transformed coordinates of an "SpatialPolygonsDataFrame" object using the projection arguments in "CRSobj", of class CRS
Because regular grids will usually not be regular after projection/datum transformation, the input object is coerced to a SpatialPointsDataFrame, and the transformation carried out on that object. A warning: “Grid warping not available, coercing to points” is given.
Because regular grids will usually not be regular after projection/datum transformation, the input object is coerced to a SpatialPointsDataFrame, and the transformation carried out on that object. A warning: “Grid warping not available, coercing to points” is given.
A non-NA logical value
The projection arguments had to be entered exactly as in the PROJ.4 documentation, in particular there cannot be any white space in +<arg>=<value> strings, and successive such strings can only be separated by blanks. Note that warnings about different projections may be issued when the PROJ.4 library extends projection arguments; examine the warning to see if the differences are real.
Also note that re-projection and/or datum transformation will usually not work for regular grids. The term used for similar operations for regular grids is warping, which involved resampling to a regular grid in the target coordinate reference system.
The methods may take an optional argument “use_ob_tran”, default FALSE, if TRUE and “+proj=ob_tran”, use General Oblique Transformation with internalised from/to projection reversal (the user oblique transforms from longlat to oblique forward rather than inverse as suggested in PROJ.4 mailing list postings); these changes are intended to meet a need pointed out by Martin Ivanov (2012-08-15). A subsequent point raised by Martin Ivanov (2017-04-28) was that use of a projected CRS with “+proj=ob_tran” led to errors, so mixing projected CRS and “+proj=ob_tran” is blocked. Transform first “+proj=ob_tran” to or from “+proj=longlat”, and then on from geographical coordinates to those desired or the reverse - see example.
If a SpatialPoints object has three dimensions, the third will also be transformed, with the metric of the third dimension assumed to be meters if the vertical units metric is not given in the projection description with +vunits= or +vto_meter= (which is 1.0 by default) https://proj.org/faq.html.
Note that WGS84 is both an ellipse and a datum, and that since 1984 there have been changes in the relative positions of continents, leading to a number of modifications. This is discussed for example in https://www.uvm.edu/giv/resources/WGS84_NAD83.pdf; there are then multiple transformations between NAD83 and WGS84 depending on the WGS84 definition used. One would expect that “+towgs84=” is a no-op for WGS84, but this only applies sometimes, and as there are now at least 30 years between now and 1984, things have shifted. It may be useful to note that “+nadgrids=@null” can help, see these threads: https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021611.html, http://lists.maptools.org/pipermail/proj/2014-August/006894.html, with thanks to Hermann Peifer for assistance.
Note that from PROJ.4 4.9.3, the definition of UTM is changed from TMERC to ETMERC; see example.
set_thin_PROJ6_warnings(TRUE)
# state
data(state)
states <- data.frame(state.x77, state.center)
states <- states[states$x > -121,]
coordinates(states) <- c("x", "y")
proj4string(states) <- CRS("+proj=longlat +ellps=clrk66")
summary(states)
#> Object of class SpatialPointsDataFrame
#> Coordinates:
#> min max
#> x -120.0680 -68.9801
#> y 27.8744 47.4231
#> Is projected: FALSE
#> proj4string : [+proj=longlat +ellps=clrk66 +no_defs]
#> Number of points: 48
#> Data attributes:
#> Population Income Illiteracy Life.Exp
#> Min. : 376 Min. :3098 Min. :0.500 Min. :67.96
#> 1st Qu.: 1188 1st Qu.:3964 1st Qu.:0.600 1st Qu.:70.13
#> Median : 2980 Median :4488 Median :0.900 Median :70.67
#> Mean : 4398 Mean :4386 Mean :1.148 Mean :70.85
#> 3rd Qu.: 5064 3rd Qu.:4768 3rd Qu.:1.450 3rd Qu.:71.88
#> Max. :21198 Max. :5348 Max. :2.800 Max. :72.96
#> Murder HS.Grad Frost Area
#> Min. : 1.400 Min. :37.80 Min. : 11.00 Min. : 1049
#> 1st Qu.: 4.275 1st Qu.:47.70 1st Qu.: 68.75 1st Qu.: 38762
#> Median : 6.850 Median :53.05 Median :114.50 Median : 54277
#> Mean : 7.321 Mean :52.64 Mean :105.65 Mean : 61749
#> 3rd Qu.:10.625 3rd Qu.:58.62 3rd Qu.:139.25 3rd Qu.: 79914
#> Max. :15.100 Max. :67.30 Max. :188.00 Max. :262134
state.ll83 <- spTransform(states, CRS("+proj=longlat +ellps=GRS80"))
summary(state.ll83)
#> Object of class SpatialPointsDataFrame
#> Coordinates:
#> min max
#> x -120.0680 -68.9801
#> y 27.8744 47.4231
#> Is projected: FALSE
#> proj4string : [+proj=longlat +ellps=GRS80 +no_defs]
#> Number of points: 48
#> Data attributes:
#> Population Income Illiteracy Life.Exp
#> Min. : 376 Min. :3098 Min. :0.500 Min. :67.96
#> 1st Qu.: 1188 1st Qu.:3964 1st Qu.:0.600 1st Qu.:70.13
#> Median : 2980 Median :4488 Median :0.900 Median :70.67
#> Mean : 4398 Mean :4386 Mean :1.148 Mean :70.85
#> 3rd Qu.: 5064 3rd Qu.:4768 3rd Qu.:1.450 3rd Qu.:71.88
#> Max. :21198 Max. :5348 Max. :2.800 Max. :72.96
#> Murder HS.Grad Frost Area
#> Min. : 1.400 Min. :37.80 Min. : 11.00 Min. : 1049
#> 1st Qu.: 4.275 1st Qu.:47.70 1st Qu.: 68.75 1st Qu.: 38762
#> Median : 6.850 Median :53.05 Median :114.50 Median : 54277
#> Mean : 7.321 Mean :52.64 Mean :105.65 Mean : 61749
#> 3rd Qu.:10.625 3rd Qu.:58.62 3rd Qu.:139.25 3rd Qu.: 79914
#> Max. :15.100 Max. :67.30 Max. :188.00 Max. :262134
state.merc <- spTransform(states, CRS=CRS("+proj=merc +ellps=GRS80"))
summary(state.merc)
#> Object of class SpatialPointsDataFrame
#> Coordinates:
#> min max
#> x -13365909 -7678830
#> y 3213175 5979931
#> Is projected: TRUE
#> proj4string :
#> [+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs]
#> Number of points: 48
#> Data attributes:
#> Population Income Illiteracy Life.Exp
#> Min. : 376 Min. :3098 Min. :0.500 Min. :67.96
#> 1st Qu.: 1188 1st Qu.:3964 1st Qu.:0.600 1st Qu.:70.13
#> Median : 2980 Median :4488 Median :0.900 Median :70.67
#> Mean : 4398 Mean :4386 Mean :1.148 Mean :70.85
#> 3rd Qu.: 5064 3rd Qu.:4768 3rd Qu.:1.450 3rd Qu.:71.88
#> Max. :21198 Max. :5348 Max. :2.800 Max. :72.96
#> Murder HS.Grad Frost Area
#> Min. : 1.400 Min. :37.80 Min. : 11.00 Min. : 1049
#> 1st Qu.: 4.275 1st Qu.:47.70 1st Qu.: 68.75 1st Qu.: 38762
#> Median : 6.850 Median :53.05 Median :114.50 Median : 54277
#> Mean : 7.321 Mean :52.64 Mean :105.65 Mean : 61749
#> 3rd Qu.:10.625 3rd Qu.:58.62 3rd Qu.:139.25 3rd Qu.: 79914
#> Max. :15.100 Max. :67.30 Max. :188.00 Max. :262134
state.merc <- spTransform(states,
CRS=CRS("+proj=merc +ellps=GRS80 +units=us-mi"))
summary(state.merc)
#> Object of class SpatialPointsDataFrame
#> Coordinates:
#> min max
#> x -8305.174 -4771.394
#> y 1996.571 3715.749
#> Is projected: TRUE
#> proj4string :
#> [+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=us-mi
#> +no_defs]
#> Number of points: 48
#> Data attributes:
#> Population Income Illiteracy Life.Exp
#> Min. : 376 Min. :3098 Min. :0.500 Min. :67.96
#> 1st Qu.: 1188 1st Qu.:3964 1st Qu.:0.600 1st Qu.:70.13
#> Median : 2980 Median :4488 Median :0.900 Median :70.67
#> Mean : 4398 Mean :4386 Mean :1.148 Mean :70.85
#> 3rd Qu.: 5064 3rd Qu.:4768 3rd Qu.:1.450 3rd Qu.:71.88
#> Max. :21198 Max. :5348 Max. :2.800 Max. :72.96
#> Murder HS.Grad Frost Area
#> Min. : 1.400 Min. :37.80 Min. : 11.00 Min. : 1049
#> 1st Qu.: 4.275 1st Qu.:47.70 1st Qu.: 68.75 1st Qu.: 38762
#> Median : 6.850 Median :53.05 Median :114.50 Median : 54277
#> Mean : 7.321 Mean :52.64 Mean :105.65 Mean : 61749
#> 3rd Qu.:10.625 3rd Qu.:58.62 3rd Qu.:139.25 3rd Qu.: 79914
#> Max. :15.100 Max. :67.30 Max. :188.00 Max. :262134
# NAD
if (FALSE) {
if (PROJis6ormore() || (!PROJis6ormore() && projNAD())) {
states <- data.frame(state.x77, state.center)
states <- states[states$x > -121,]
coordinates(states) <- c("x", "y")
proj4string(states) <- CRS("+init=epsg:4267")
print(summary(states))
state.ll83 <- spTransform(states, CRS("+init=epsg:4269"))
print(summary(state.ll83))
state.kansasSlcc <- spTransform(states, CRS=CRS("+init=epsg:26978"))
print(summary(state.kansasSlcc))
SFpoint_NAD83 <- SpatialPoints(matrix(c(-103.869667, 44.461676), nrow=1),
proj4string=CRS("+init=epsg:4269"))
SFpoint_NAD27 <- spTransform(SFpoint_NAD83, CRS("+init=epsg:4267"))
print(all.equal(coordinates(SFpoint_NAD83), coordinates(SFpoint_NAD27)))
print(coordinates(SFpoint_NAD27), digits=12)
print(coordinates(SFpoint_NAD83), digits=12)
}
}
data(meuse)
coordinates(meuse) <- c("x", "y")
proj4string(meuse) <- CRS("+init=epsg:28992")
# see http://trac.osgeo.org/gdal/ticket/1987
summary(meuse)
#> Object of class SpatialPointsDataFrame
#> Coordinates:
#> min max
#> x 178605 181390
#> y 329714 333611
#> Is projected: TRUE
#> proj4string :
#> [+proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889
#> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs]
#> Number of points: 155
#> Data attributes:
#> cadmium copper lead zinc
#> Min. : 0.200 Min. : 14.00 Min. : 37.0 Min. : 113.0
#> 1st Qu.: 0.800 1st Qu.: 23.00 1st Qu.: 72.5 1st Qu.: 198.0
#> Median : 2.100 Median : 31.00 Median :123.0 Median : 326.0
#> Mean : 3.246 Mean : 40.32 Mean :153.4 Mean : 469.7
#> 3rd Qu.: 3.850 3rd Qu.: 49.50 3rd Qu.:207.0 3rd Qu.: 674.5
#> Max. :18.100 Max. :128.00 Max. :654.0 Max. :1839.0
#>
#> elev dist om ffreq soil lime
#> Min. : 5.180 Min. :0.00000 Min. : 1.000 1:84 1:97 0:111
#> 1st Qu.: 7.546 1st Qu.:0.07569 1st Qu.: 5.300 2:48 2:46 1: 44
#> Median : 8.180 Median :0.21184 Median : 6.900 3:23 3:12
#> Mean : 8.165 Mean :0.24002 Mean : 7.478
#> 3rd Qu.: 8.955 3rd Qu.:0.36407 3rd Qu.: 9.000
#> Max. :10.520 Max. :0.88039 Max. :17.000
#> NA's :2
#> landuse dist.m
#> W :50 Min. : 10.0
#> Ah :39 1st Qu.: 80.0
#> Am :22 Median : 270.0
#> Fw :10 Mean : 290.3
#> Ab : 8 3rd Qu.: 450.0
#> (Other):25 Max. :1000.0
#> NA's : 1
meuse.utm <- spTransform(meuse, CRS("+proj=utm +zone=32 +datum=WGS84"))
#> Warning: PROJ support is provided by the sf and terra packages among others
summary(meuse.utm)
#> Object of class SpatialPointsDataFrame
#> Coordinates:
#> min max
#> x 269909.9 272832.1
#> y 5650085.3 5653888.9
#> Is projected: TRUE
#> proj4string :
#> [+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs]
#> Number of points: 155
#> Data attributes:
#> cadmium copper lead zinc
#> Min. : 0.200 Min. : 14.00 Min. : 37.0 Min. : 113.0
#> 1st Qu.: 0.800 1st Qu.: 23.00 1st Qu.: 72.5 1st Qu.: 198.0
#> Median : 2.100 Median : 31.00 Median :123.0 Median : 326.0
#> Mean : 3.246 Mean : 40.32 Mean :153.4 Mean : 469.7
#> 3rd Qu.: 3.850 3rd Qu.: 49.50 3rd Qu.:207.0 3rd Qu.: 674.5
#> Max. :18.100 Max. :128.00 Max. :654.0 Max. :1839.0
#>
#> elev dist om ffreq soil lime
#> Min. : 5.180 Min. :0.00000 Min. : 1.000 1:84 1:97 0:111
#> 1st Qu.: 7.546 1st Qu.:0.07569 1st Qu.: 5.300 2:48 2:46 1: 44
#> Median : 8.180 Median :0.21184 Median : 6.900 3:23 3:12
#> Mean : 8.165 Mean :0.24002 Mean : 7.478
#> 3rd Qu.: 8.955 3rd Qu.:0.36407 3rd Qu.: 9.000
#> Max. :10.520 Max. :0.88039 Max. :17.000
#> NA's :2
#> landuse dist.m
#> W :50 Min. : 10.0
#> Ah :39 1st Qu.: 80.0
#> Am :22 Median : 270.0
#> Fw :10 Mean : 290.3
#> Ab : 8 3rd Qu.: 450.0
#> (Other):25 Max. :1000.0
#> NA's : 1
cbind(coordinates(meuse), coordinates(meuse.utm))
#> x y x y
#> 1 181072 333611 272531.7 5653889
#> 2 181025 333558 272482.1 5653838
#> 3 181165 333537 272620.9 5653810
#> 4 181298 333484 272751.2 5653751
#> 5 181307 333330 272752.7 5653597
#> 6 181390 333260 272832.1 5653523
#> 7 181165 333370 272612.8 5653644
#> 8 181027 333363 272474.5 5653643
#> 9 181060 333231 272501.0 5653510
#> 10 181232 333168 272669.8 5653438
#> 11 181191 333115 272626.2 5653388
#> 12 181032 333031 272463.3 5653311
#> 13 180874 333339 272320.5 5653627
#> 14 180969 333252 272411.2 5653535
#> 15 181011 333161 272448.7 5653442
#> 16 180830 333246 272272.0 5653536
#> 17 180763 333104 272198.1 5653398
#> 18 180694 332972 272122.7 5653269
#> 19 180625 332847 272047.6 5653148
#> 20 180555 332707 271970.8 5653011
#> 21 180642 332708 272057.8 5653008
#> 22 180704 332717 272120.2 5653014
#> 23 180704 332664 272117.6 5652961
#> 24 181153 332925 272578.9 5653200
#> 25 181147 332823 272567.9 5653098
#> 26 181167 332778 272585.7 5653052
#> 27 181008 332777 272426.8 5653059
#> 28 180973 332687 272387.4 5652971
#> 29 180916 332753 272333.7 5653039
#> 30 181352 332946 272778.8 5653211
#> 31 181133 332570 272541.5 5652846
#> 32 180878 332489 272282.8 5652777
#> 33 180829 332450 272232.0 5652741
#> 34 180954 332399 272354.3 5652684
#> 35 180956 332318 272352.4 5652603
#> 37 180710 332330 272107.2 5652627
#> 38 180632 332445 272034.9 5652746
#> 39 180530 332538 271937.6 5652843
#> 40 180478 332578 271887.6 5652886
#> 41 180383 332476 271787.7 5652789
#> 42 180494 332330 271891.4 5652637
#> 43 180561 332193 271951.6 5652497
#> 44 180451 332175 271840.8 5652485
#> 45 180410 332031 271792.8 5652343
#> 46 180355 332299 271751.0 5652613
#> 47 180292 332157 271681.1 5652475
#> 48 180283 332014 271665.1 5652332
#> 49 180282 331861 271656.6 5652179
#> 50 180270 331707 271637.1 5652026
#> 51 180199 331591 271560.5 5651914
#> 52 180135 331552 271494.6 5651878
#> 53 180237 332351 271635.7 5652671
#> 54 180103 332297 271499.2 5652624
#> 55 179973 332255 271367.2 5652588
#> 56 179826 332217 271218.5 5652557
#> 57 179687 332161 271076.9 5652508
#> 58 179792 332035 271175.6 5652377
#> 59 179902 332113 271289.3 5652450
#> 60 180100 332213 271492.1 5652540
#> 61 179604 332059 270989.0 5652410
#> 62 179526 331936 270905.0 5652291
#> 63 179495 331770 270865.9 5652127
#> 64 179489 331633 270853.2 5651990
#> 65 179414 331494 270771.5 5651855
#> 66 179334 331366 270685.3 5651731
#> 67 179255 331264 270601.3 5651633
#> 69 179470 331125 270809.3 5651484
#> 75 179692 330933 271021.7 5651281
#> 76 179852 330801 271175.1 5651141
#> 79 179140 330955 270471.3 5651330
#> 80 179128 330867 270455.0 5651243
#> 81 179065 330864 270391.9 5651243
#> 82 179007 330727 270327.3 5651109
#> 83 179110 330758 270431.7 5651135
#> 84 179032 330645 270348.2 5651026
#> 85 179095 330636 270410.7 5651014
#> 86 179058 330510 270367.6 5650890
#> 87 178810 330666 270127.5 5651058
#> 88 178912 330779 270234.9 5651166
#> 89 178981 330924 270310.9 5651307
#> 90 179076 331005 270409.8 5651383
#> 123 180151 330353 271451.8 5650679
#> 160 179211 331175 270553.0 5651547
#> 163 181118 333214 272558.1 5653490
#> 70 179474 331304 270822.1 5651662
#> 71 179559 331423 270912.8 5651777
#> 91 179022 330873 270349.4 5651254
#> 92 178953 330742 270274.0 5651127
#> 93 178875 330516 270185.0 5650905
#> 94 178803 330349 270104.9 5650741
#> 95 179029 330394 270332.9 5650775
#> 96 178605 330406 269909.9 5650808
#> 97 178701 330557 270013.2 5650954
#> 98 179547 330245 270843.1 5650601
#> 99 179301 330179 270594.1 5650547
#> 100 179405 330567 270717.0 5650930
#> 101 179462 330766 270783.7 5651126
#> 102 179293 330797 270616.4 5651165
#> 103 179180 330710 270499.3 5651083
#> 104 179206 330398 270509.9 5650771
#> 105 179618 330458 270924.5 5650810
#> 106 179782 330540 271092.3 5650884
#> 108 179980 330773 271301.6 5651107
#> 109 180067 331185 271408.7 5651515
#> 110 180162 331387 271513.5 5651712
#> 111 180451 331473 271806.4 5651783
#> 112 180328 331158 271668.1 5651475
#> 113 180276 330963 271606.6 5651283
#> 114 180114 330803 271436.9 5651131
#> 115 179881 330912 271209.5 5651251
#> 116 179774 330921 271103.0 5651265
#> 117 179657 331150 270997.4 5651500
#> 118 179731 331245 271075.9 5651591
#> 119 179717 331441 271071.6 5651787
#> 120 179446 331422 270799.9 5651782
#> 121 179524 331565 270884.8 5651921
#> 122 179644 331730 271012.8 5652080
#> 124 180321 330366 271622.3 5650684
#> 125 180162 331837 271535.6 5652161
#> 126 180029 331720 271396.9 5652051
#> 127 179797 331919 271174.9 5652261
#> 128 179642 331955 271021.9 5652305
#> 129 179849 332142 271237.8 5652481
#> 130 180265 332297 271661.0 5652616
#> 131 180107 332101 271493.6 5652428
#> 132 180462 331947 271840.7 5652256
#> 133 180478 331822 271850.5 5652131
#> 134 180347 331700 271713.7 5652015
#> 135 180862 333116 272297.6 5653405
#> 136 180700 332882 272124.3 5653179
#> 161 180201 331160 271541.3 5651483
#> 162 180173 331923 271550.8 5652247
#> 137 180923 332874 272346.7 5653160
#> 138 180467 331694 271833.2 5652003
#> 140 179917 331325 271265.7 5651662
#> 141 179822 331242 271166.7 5651583
#> 142 179991 331069 271327.1 5651402
#> 143 179120 330578 270432.8 5650955
#> 144 179034 330561 270346.1 5650942
#> 145 179085 330433 270390.8 5650811
#> 146 179236 330046 270522.6 5650417
#> 147 179456 330072 270743.7 5650433
#> 148 179550 329940 270831.1 5650296
#> 149 179445 329807 270719.7 5650168
#> 150 179337 329870 270614.9 5650237
#> 151 179245 329714 270515.4 5650085
#> 152 179024 329733 270295.5 5650115
#> 153 178786 329822 270062.1 5650216
#> 154 179135 329890 270414.1 5650267
#> 155 179030 330082 270318.6 5650463
#> 156 179184 330182 270477.4 5650556
#> 157 179085 330292 270383.9 5650671
#> 158 178875 330311 270175.0 5650700
#> 159 179466 330381 270768.8 5650741
#> 164 180627 330190 271919.3 5650493
if (FALSE) {
# Kiritimati
kiritimati_primary_roads <- readOGR(system.file("vectors",
package = "rgdal")[1], "kiritimati_primary_roads")
kiritimati_primary_roads_ll <- spTransform(kiritimati_primary_roads,
CRS("+proj=longlat +datum=WGS84"))
opar <- par(mfrow=c(1,2))
plot(kiritimati_primary_roads, axes=TRUE)
plot(kiritimati_primary_roads_ll, axes=TRUE, las=1)
par(opar)
# scot_BNG
scot_BNG <- readOGR(system.file("vectors", package = "rgdal")[1],
"scot_BNG")
scot_LL <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84"))
plot(scot_LL, axes=TRUE)
grdtxt_LL <- gridat(scot_LL)
grd_LL <- gridlines(scot_LL, ndiscr=100)
summary(grd_LL)
target <- CRS(proj4string(scot_BNG))
grd_BNG <- spTransform(grd_LL, target)
grdtxt_BNG <- spTransform(grdtxt_LL, target)
opar <- par(mfrow=c(1,2))
plot(scot_BNG, axes=TRUE, las=1)
plot(grd_BNG, add=TRUE, lty=2)
text(coordinates(grdtxt_BNG),
labels=parse(text=as.character(grdtxt_BNG$labels)))
par(opar)
}
# broke_proj
broke_proj <- FALSE
# https://github.com/OSGeo/PROJ/issues/1525
pv <- .Call("PROJ4VersionInfo", PACKAGE="rgdal")[[2]]
if (pv >= 600 && pv < 620) broke_proj <- TRUE
if (!broke_proj) {
crds <- matrix(data=c(9.05, 48.52), ncol=2)
spPoint <- SpatialPoints(coords=crds,
proj4string=CRS("+proj=longlat +ellps=sphere +no_defs"))
ob_tran_def <- paste("+proj=ob_tran +o_proj=longlat",
"+o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs")
tg <- CRS(ob_tran_def)
# proj4string not propagated in GDAL 3.0.0
a <- spTransform(spPoint, tg, use_ob_tran=TRUE)
a
}
#> SpatialPoints:
#> coords.x1 coords.x2
#> [1,] -5.917698 -1.87195
#> Coordinate Reference System (CRS) arguments: +proj=ob_tran
#> +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere
#> +no_defs
#should be (-5.917698, -1.87195)
if (!broke_proj) {
spTransform(a, CRS("+proj=longlat +ellps=sphere +no_defs"),
use_ob_tran=TRUE)
}
#> SpatialPoints:
#> coords.x1 coords.x2
#> [1,] 9.05 48.52
#> Coordinate Reference System (CRS) arguments: +proj=longlat
#> +ellps=sphere +no_defs
if (!broke_proj) {
try(spTransform(a, CRS(paste("+proj=tmerc +lat_0=0 +lon_0=9 +k=1",
"+x_0=3500000 +y_0=0 +ellps=bessel +units=m +no_defs")),
use_ob_tran=TRUE))
}
#> Error in spTransform(a, CRS(paste("+proj=tmerc +lat_0=0 +lon_0=9 +k=1", :
#> Use ob_tran to or from unprojected objects only
if (!broke_proj) {
spTransform(spPoint, CRS(paste("+proj=tmerc +lat_0=0 +lon_0=9 +k=1",
"+x_0=3500000 +y_0=0 +ellps=bessel +units=m +no_defs")))
}
#> SpatialPoints:
#> coords.x1 coords.x2
#> [1,] 3503693 5375702
#> Coordinate Reference System (CRS) arguments: +proj=tmerc +lat_0=0
#> +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +units=m +no_defs
if (!broke_proj) {
spTransform(spTransform(a, CRS("+proj=longlat +ellps=sphere +no_defs"),
use_ob_tran=TRUE), CRS(paste("+proj=tmerc +lat_0=0 +lon_0=9 +k=1",
"+x_0=3500000 +y_0=0 +ellps=bessel +units=m +no_defs")))
}
#> SpatialPoints:
#> coords.x1 coords.x2
#> [1,] 3503693 5375702
#> Coordinate Reference System (CRS) arguments: +proj=tmerc +lat_0=0
#> +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +units=m +no_defs
crds1 <- matrix(data=c(7, 51, 8, 52, 9, 52, 10, 51, 7, 51), ncol=2,
byrow=TRUE, dimnames=list(NULL, c("lon", "lat")));
crds2 <- matrix(data=c(8, 48, 9, 49, 11, 49, 9, 48, 8, 48), ncol=2,
byrow=TRUE, dimnames=list(NULL, c("lon", "lat")));
crds3 <- matrix(data=c(6, 47, 6, 55, 15, 55, 15, 47, 6, 47), ncol=2,
byrow=TRUE, dimnames=list(NULL, c("lon", "lat")));
spLines <- SpatialLines(list(Lines(list(Line(crds1), Line(crds2),
Line(crds3)), ID="a")));
slot(spLines, "proj4string") <- CRS("+proj=longlat +ellps=sphere +no_defs");
bbox(spLines);
#> min max
#> x 6 15
#> y 47 55
if (!broke_proj) {
spLines_tr <- spTransform(spLines, tg, use_ob_tran=TRUE);
bbox(spLines_tr)
}
#> min max
#> x -8.163685 -1.725025
#> y -3.708442 4.807871
if (!broke_proj) {
bbox(spTransform(spLines_tr, CRS("+proj=longlat +ellps=sphere"),
use_ob_tran=TRUE))
}
#> min max
#> x 6 15
#> y 47 55
if (!broke_proj) {
spPolygons <- SpatialPolygons(list(Polygons(list(Polygon(crds1),
Polygon(crds2), Polygon(crds3)), ID="a")));
slot(spPolygons, "proj4string") <- CRS("+proj=longlat +ellps=sphere +no_defs");
bbox(spPolygons);
}
#> min max
#> x 6 15
#> y 47 55
if (!broke_proj) {
spPolygons_tr <- spTransform(spPolygons, tg, use_ob_tran=TRUE);
bbox(spPolygons_tr)
}
#> min max
#> x -8.163685 -1.725025
#> y -3.708442 4.807871
if (!broke_proj) {
bbox(spTransform(spPolygons_tr, CRS("+proj=longlat +ellps=sphere"),
use_ob_tran=TRUE))
}
#> min max
#> x 6 15
#> y 47 55
# added after posting by Martin Ivanov
if (FALSE) {
data(nor2k)
summary(nor2k)
nor2kNGO <- spTransform(nor2k, CRS("+init=epsg:4273"))
summary(nor2kNGO)
all.equal(coordinates(nor2k)[,3], coordinates(nor2kNGO)[,3])
# added after posting by Don MacQueen
crds <- cbind(c(-121.524764291826, -121.523480804667), c(37.6600366036405, 37.6543604613483))
ref <- cbind(c(1703671.30566227, 1704020.20113366), c(424014.398045834, 421943.708664294))
crs.step1.cf <- CRS(paste("+proj=lcc +lat_1=38.43333333333333",
"+lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5",
"+x_0=2000000.0 +y_0=500000.0 +ellps=GRS80 +units=us-ft +no_defs",
"+towgs84=-0.991,1.9072,0.5129,0.025789908,0.0096501,0.0116599,0.0"))
locs.step1.cf <- spTransform(SpatialPoints(crds,
proj4string=CRS("+proj=longlat +datum=WGS84")), crs.step1.cf)
suppressWarnings(proj4string(locs.step1.cf) <- CRS(paste("+proj=lcc",
"+lat_1=38.43333333333333 +lat_2=37.06666666666667 +lat_0=36.5",
"+lon_0=-120.5 +x_0=2000000.0 +y_0=500000.0 +ellps=GRS80 +units=us-ft",
"+no_defs +nadgrids=@null")))
locs.step2.cfb <- spTransform(locs.step1.cf, CRS("+init=epsg:26743"))
coordinates(locs.step2.cfb) - ref
all.equal(unname(coordinates(locs.step2.cfb)), ref)
}
if (FALSE) {
# new_proj_and_gdal()
run <- new_proj_and_gdal()
if (run) {
# Test for UTM == TMERC (<= 4.9.2) or UTM == ETMERC (> 4.9.2)
nhh <- SpatialPointsDataFrame(matrix(c(5.304234, 60.422311), ncol=2),
proj4string=CRS(SRS_string="OGC:CRS84"), data=data.frame(office="RSB"))
nhh_utm_32N_P4 <- spTransform(nhh, CRS("+init=epsg:3044"))
nhh_tmerc_P4 <- spTransform(nhh, CRS(paste("+proj=tmerc +k=0.9996",
"+lon_0=9 +x_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")))
nhh_etmerc_P4 <- spTransform(nhh, CRS(paste("+proj=etmerc +k=0.9996",
"+lon_0=9 +x_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")))
all.equal(coordinates(nhh_utm_32N_P4), coordinates(nhh_tmerc_P4),
tolerance=1e-9, scale=1)
# UTM == TMERC: PROJ4 <=4.9.2
all.equal(coordinates(nhh_utm_32N_P4), coordinates(nhh_etmerc_P4),
tolerance=1e-9, scale=1)
# UTM == ETMERC: PROJ4 > 4.9.2
unis <- SpatialPointsDataFrame(matrix(c(15.653453, 78.222504), ncol=2),
proj4string=CRS(SRS_string="OGC:CRS84"), data=data.frame(office="UNIS"))
unis_utm_33N_P4 <- spTransform(unis, CRS("+init=epsg:3045"))
unis_tmerc_P4 <- spTransform(unis, CRS(paste("+proj=tmerc +k=0.9996 +lon_0=15",
"+x_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")))
unis_etmerc_P4 <- spTransform(unis, CRS(paste("+proj=etmerc +k=0.9996",
"+lon_0=15 +x_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")))
all.equal(coordinates(unis_utm_33N_P4), coordinates(unis_tmerc_P4),
tolerance=1e-9, scale=1)
# UTM == TMERC: PROJ4 <=4.9.2
all.equal(coordinates(unis_utm_33N_P4), coordinates(unis_etmerc_P4),
tolerance=1e-9, scale=1)
# UTM == ETMERC: PROJ4 > 4.9.2
}
}