Transform Coordinates to a Different Projected CRS
If you directly compare data sets with different projected coordinate reference systems (CRSs), then the results are inaccurate. Therefore, before comparing data sets, first verify that the CRSs are the same. If different projected CRSs have the same underlying geographic CRS, then you can transform the coordinates from one projected CRS to the other. Once the data sets are referenced to the same projected CRS, you can compare them.
将projectedx-ycoordinates to a different projected CRS, first unproject thex-ycoordinates to latitude-longitude coordinates by using theprojinv
function. Then, project the latitude-longitude coordinates tox-ycoordinates in a different projected CRS by using theprojfwd
function.
For example, import a shapefile containing thex- andy-coordinates of roads in Boston. Also import information about the shapefile as a structure. Find the projected CRS for the coordinates by accessing theCoordinateReferenceSystem
field of that structure.
s = shaperead('boston_roads.shp'); x1 = [s.X]; y1 = [s.Y]; info = shapeinfo('boston_roads.shp'); p1 = info.CoordinateReferenceSystem;
Unproject thex-ycoordinates and return latitude-longitude coordinates.
[lat,lon] = projinv(p1,x1,y1);
Select a new projected CRS for the target projection. For this example, create aprojcrs
object for UTM zone 19N. Verify that both projected CRSs have the same geographic CRS. If the geographic CRSs are different, then the projected coordinates may be inaccurate. You can find the geographic CRS by querying theGeographicCRS
property of theprojcrs
object.
p2 = projcrs(26919); p2.GeographicCRS.Name
ans = "NAD83"
p1.GeographicCRS.Name
ans = "NAD83"
Project the latitude-longitude coordinates tox-ycoordinates by specifying theprojcrs
object you created.
[x2,y2] = projfwd(p2,lat,lon);
Compare the originalx-ycoordinates with the newx-ycoordinates by displaying them. Add labels and a title to each figure.
figure mapshow(x1,y1) xlabel('x (meters)') ylabel('y (meters)') title(p1.Name)
figure mapshow(x2,y2) xlabel('x (meters)') ylabel('y (meters)') title(p2.Name)
The visualizations are similar, but the coordinates displayed along the axis rulers correspond to different projected CRSs.
See Also
projinv
|projfwd
|geocrs
|projcrs
|shaperead