How do I change the projection of raster data in MATLAB?

Is it possible to change the projection of georeferenced raster data in Matlab? I am interested in both reprojecting the my spatial data (geotiffs) between different projected coordinate systems and "unprojecting" it to geographic coordinates. Basically, I'm looking to do the equivalent "project raster" in ArcGIS. I can project and un-project shapefiles without issue using minvtran and mfwdtran, but I can't find an equivalent function for rasters.

Thanks!

Answer 1

[Chad Greene]

Halushka,

The mfwdtran and minvtran functions will work. If you have some x and y vectors associated with your geotiff, grid the coordinates like this:

[X,Y] = meshgrid(x,y);

then use minvtran to get geocoordinates of each pixel by

[lat,lon] = minvtran(X,Y);

Then you'll be able to use pcolorm(lat,lon,Z) or transform lat and lon into a different projection using mfwdtran.

[Halushka]

Thanks for answering Chad. Sorry for the delay in my response, but I wanted to take the time to sit down with your answer and try to sort my issue out on my own. I finally had the time, but unfortunately I am still confused:

As an example, let's say I have gridded climate data in geotiff format in a lambert conformal conic projection. Let's say I also have shapefiles bounding several regions of interest and I would like to analyze the climate data by region. The first thing I need to do is reproject the data to an equal-area projection. I know I can do this:

 ras = 'MyRaster_Lambert.tif';

% Load geotiff
[geotiff,R] = geotiffread(ras);
[info] = geotiffinfo(ras);

% Get pixel centers
[x,y] = pixcenters(R,size(geotiff));
[X,Y] = meshgrid(x,y);

% Get geocoordinates
mstruct = geotiff2mstruct(info);
[lat,lon] = minvtran(mstruct,X,Y);

From your answer, I figured out I can further do this:

P = 'eqacylin'; % an equal-area projection
E = referenceEllipsoid('wgs84');
axesm('MapProjection',P,'Geoid',E);
h = pcolorm(lat,lon,geotiff);

I now have a my geotiff displayed as in an equal area projection. BUT...

How do I get at the actual data? The next step would be for me to mask my raster data with the shapefiles delineating my regions of interest (this I have been able to do using inpolygon as long as my input data and shapefiles were already in the same projection). To do that I need my climate data in a matrix as it is in geotiff, and I need to construct a new R struct to get the new pixel center coordinates. This is where I am stuck. Any thoughts would be appreciated, I'm (obviously) new to the Mapping Toolbox and I am not finding it very intuitive.

[guruprasad]

Hey, were you able to figure out how to get the actual reprojected data? I am stuck with a similar issue.

Cheers, Guru

[Halushka]

No I wasn't.

Answer 2

[Matthew Cooper]

I recommend using gdal, which can be called from Matlab using the system command or the shell escape character. Here is an example that pretends you want to reproject your data from North America Lambert Conformal Conic projection to Polar Stereographic North projection:

ras_in  = 'MyRaster_Lambert.tif';
ras_out = 'MyRaster_Stereo.tif';
func    = '/path/to/gdalwarp ';
opts    = '-s_srs ESRI:102009 -t_srs EPSG:3413 -r bilinear -of GTiff -co "TFW=YES" ';
command = [func opts ras_in ' ' ras_out];
status  = system(command);

You will have to replace the path to gdalwarp with the path where it exists on your system. Mine is /usr/local/bin/gdalwarp. You can customize this by defining path variables or adding the path to your matlab search path. In the 'opts' variable, I included some common gdalwarp options that might be helpful, but check out the gdal documention to understand what is available.

The same command using the shell escape character would look like this:

!/path/to/gdalwarp -s_srs ESRI:102009 -t_srs EPSG:3413 -r bilinear -of GTiff -co "TFW=YES" MyRaster_Lambert.tif MyRaster_Stereo.tif

You can run the command above directly in the matlab editor/terminal.

Performing the reprojection with Matlab is possible but quite difficult and involves careful thought, way beyond the scope of this answer. In its simplest form, you need to interpolate the original data onto the new coordinates, but you need to consider what method of interpolation is appropriate and how the projection affects the mapping between the original data and the 'projected' data.

One last point. Matlab will not allow you to change the R.RasterInterpretation property between 'Cells' and 'Postings'. This has important implications for all of the other functions in the mapping toolbox. Make sure you understand which interpretation is correct. If you need to modify the interpretation before reading the data into matlab, you can use something like this:

!/path/to/gdal_translate -mo AREA_OR_POINT=AREA MyRaster_Stereo.tif MyRaster_Stereo_asCells.tif