Search code examples
matlabcoordinatesrastermap-projections

Projecting lat-lon coordinates to a raster to obtain pixel values at those coordinates

We have a georeferenced raster file and a list of latitude-longitude coordinates that we want to plot on top of that raster file to obtain the raster values at those coordinates.

Here it is the list of coordinates:

mapshow(lon,lat,'DisplayType','point')  %x,y

enter image description here

Of course, longitude is between -180 and +180 and latitude between -90 and +90.

And here it is the raster file:

proj = geotiffinfo('global2001.tif');
[raster_tiff,cmap_tiff,reference_tiff] = geotiffread('global2001.tif');
figure
mapshow(raster_tiff,cmap_tiff,reference_tiff)

enter image description here

proj =           FileModDate: '20-nov-2018 14:15:52'
             FileSize: 121625752
               Format: 'tif'
        FormatVersion: []
               Height: 40032
                Width: 80062
             BitDepth: 8
            ColorType: 'indexed'
            ModelType: 'ModelTypeProjected'
                  PCS: ''
           Projection: ''
               MapSys: ''
                 Zone: []
         CTProjection: 'CT_Sinusoidal'
             ProjParm: [7×1 double]
           ProjParmId: {7×1 cell}
                  GCS: 'WGS 84'
                Datum: 'World Geodetic System 1984'
            Ellipsoid: 'WGS 84'
            SemiMajor: 6378137
            SemiMinor: 6.3568e+06
                   PM: 'Greenwich'
    PMLongToGreenwich: 0
            UOMLength: 'metre'
    UOMLengthInMeters: 1
             UOMAngle: 'degree'
    UOMAngleInDegrees: 1
            TiePoints: [1×1 struct]
           PixelScale: [3×1 double]
           SpatialRef: [1×1 map.rasterref.MapCellsReference]
            RefMatrix: [3×2 double]
          BoundingBox: [2×2 double]
         CornerCoords: [1×1 struct]
         GeoTIFFCodes: [1×1 struct]
          GeoTIFFTags: [1×1 struct]

And now we project the coordinates using the same projection as the raster file:

mstruct = geotiff2mstruct(proj);
% get current axis 
axesm('MapProjection',mstruct.mapprojection,'Frame','on')
h=get(gcf,'CurrentAxes');
assert(ismap(h)==1,'Current axes must be map axes.')
mstruct=gcm;
[x,y] = mfwdtran(mstruct,lat,lon,h,'none');

figure
mapshow(x,y,'DisplayType','point')  %x,y

Resulting in this figure enter image description here

The problem is that the projected coordinates go from -3 to +3 and from -1 to +1 and the axis in the raster file goes from -2 to +2 but to the power of 7, so if we plot those points on top of that raster we see everything as one single point somewhere in the Atlantic.

Eventually, we want to use the function latlon2pix in order to obtain the pixel values at each of the coordinate points, but first, we need to be capable of puting the two things together. Any ideas?

The geoshow function doesn't work. We have enough RAM...

Solution

  • I suggest this:

    % assigning reference for tif file    
proj = geotiffinfo('global2001.tif');
% reading the tif into raster_tiff, storing the colormap, save the meta data 
[raster_tiff,cmap_tiff,reference_tiff] = geotiffread('global2001.tif');

    The projection can be done using projfwd and the projection matrix stored in variable proj. Data points will then be expressed in (here sinusoidal) map coordinates of the raster file.

    % using the projection matrix (proj.RefMatrix) and computing the raster coordinates in meters
[x,y] = projfwd(proj,lat,lon);

    Using mapshow one can plot the two data sets on top. 

    mapshow(raster_tiff,cmap_tiff,reference_tiff);
mapshow(x,y,'DisplayType', 'point', 'Color', 'm',...
            'MarkerEdgeColor', 'auto');