How to change coordinate projection system of a shapfile

Using the MATLAB Mapping Toolbox, I am working with raster (Landsat GeoTIFF) and vector data (externally provided shapefiles) that have different coordinate projection systems. The GeoTIFFs are in WGS84_UTM_zone_50S (PCS = 32750) while the shapefiles are in AGD84_AMG_zone_50 (PCS = 20350). Hence, when overlaying the shapefile polygons on the raster image they are out by about 200m.

One way to fix this could be to rewrite the GeoTIFFs and alter their PCS to AGD84_AMG_zone_50 and work in that PCS, but I'd prefer to alter the PCS of the shapefile. I've read about mapstruct and mstruct but the fog hasn't cleared - if there's a way of using these (or other functions) to achieve my goal then I'd really appreciate a simple explanation.

ARCGIS can project data sources with differing PCS 'on-the-fly' but I guess the Mapping Toolbox can't do this?

Kind regards, Craig.

Answer 1

[Craig]

Thanks for the ideas Chad,

I tried the two approaches above and failed but I did learn a lot. The problem is that the shapefile is in map units and, although it has a projected coordinate system (PCS), MATLAB doesn't read the optional '.prj' file associated with the shapefile and doesn't build a projection structure for the shapefile data. So:

1. Using geoshow isn't possible because the units are mapping units.

2. "'UseGeoCoords', 'true'" doesn't convert mapping units to geographic units - it just tells the function that the data are already in geographic units when reading them in.

3. The projinv and related functions won't work unless the user manually defines a projection structure for the data that MATLAB will recognise. This is possible but a lot of work.

It should be possible to alter the GeoTIFF's PCS using projinv and projfwd because it is read in with a projection structure. But I don't really want to alter that. My solution at the moment is to alter the shapefile's PCS in ARCGIS and save as a new Feature Class using the Project Tool.

It irks me having to leave the MATLAB environment and use ARCGIS for this simple transformation but there it is. Hopefully a future version of the Mapping Toolbox will read in the '.prj' file info and convert to a projection structure for use with the toolboxes functions.

Thanks for helping me to elucidate the problem if not solve it.

Answer 2

[Chad Greene]

I've never mixed coordinate systems, but it doesn't seem like it'd be a problem to plot one dataset using geoshow or mapshow with its proper spatial referencing object R, then overlay some other data set that may or may not have the same spatial referencing object.

Another thought: I have not worked with shapefiles, but shaperead has a UseGeoCoords options, which returns data in lat/lon format, independent of any projection information. This seems like it should work:

[A, R] = geotiffread('someLandsatImage.TIF'); % loads image
worldmap('world') % initializes map 
geoshow(A,R) % projects Landsat image
S = shaperead('worldcities.shp', 'UseGeoCoords', true); % imports lat/lon data
textm(S.Lat,S.Lon,S.Name) % puts lat/lon data on map

[Alexander Audet]

As the Asker mentioned,

'UseGeoCoords', true

does not convert your coordinates to GeoCoordinates, it just tells Matlab (mistakenly in this case) that they are GeoCoordinates. This answer is misleading and should be changed.

Answer 3

[Chad Greene]

Can you unproject the polygon data with projinv, then use projfwd to project into the Landsat projection?

Answer 4

[Chad Greene]

Do any of these tools help?

[Craig]

I think it would be possible if I had, or were willing to obtain, a deep enough understanding of geodetic transformations. But it would take a lot of time and what I love about MATLAB is how quick and easy it makes most things - but this isn't one of them. I can alter the PCS in ARCGIS in about five minutes so I'll stick with that for now.

But thanks for bringing that toolbox to my notice.

Answer 5

[Chad Greene]

There's no magic to coordinate transformations. If you have to perform the transformation more than once, it might be worth it to create some simple conversion functions. And if they help you, they'll probably help someone else, so share them on the File Exchange :) I find that when I share my functions on FEX, my code ends up being more efficient and better documented than when I write one-time-use scripts. Also, the process of documenting code creates a much deeper understanding of the formulas and your own code. Of course you can do the ArcGIS conversions, but then you'd have to use ArcGIS. Personally, I'd rather spend an hour writing a Matlab function than spend a minute using ArcGIS.