matlab与geotiff影像的函数，用法介绍_pix2map

---

前言

matlab给影像加地理坐标信息，一直没有很系统的明白，现在学习并记录吧

---

一、matlab和geotiff相关函数的学习

(1) pix2latlon() 和latlon2pix()

pix2latlon()和latlon2pix()是对经纬度坐标系的影像，实现**像素坐标**和**经纬度坐标**的转换

%% (1)pix2latlon( )        [lat,lon]=pix2latlon(R,row,col) 
%     latlon2pix()           [row,col]=pix2latlon(R,lat,lon) 
%自己的理解： 应该是对地理坐标系的影像，进行像素坐标系和地理坐标系的转换。
% R是 geotiff影像的info.RefMatrix(3行2列)，或者是一个 a geographic raster reference object.
img_name='E:\2022_5_2_GF_pa1_pa2_lens_rect\pa1_012_006\mapproject\pm1_006_6087_2_12116_sub_ba_rect.tif';
row=100;col=100;
[A,refmat,bbox] = geotiffread(img_name);% geotiffread的第二个参数就是geotif 的info.RefMatrix
info = geotiffinfo(img_name);
R=info.RefMatrix
[lat,lon]=pix2latlon(R,row,col);
% [lat,lon]=pix2latlon(refmat,row,col);
% R的读取，可以用过geotiffread，也可以通过geotiffinfo，info.RefMatrix

1
2
3
4
5
6
7
8
9
10
11
12
13

(2) pix2map() 和 map2pix()

（1）pix2map()和map2pix()是对投影坐标系的影像，实现**投影坐标**<--->**像素坐标**的转换
（2）当然，也可以实现**像素坐标**和**经纬度坐标**的转换，需要借用projfwd（）和projinv（）函数

%% (2)pix2map()           [x,y] = pix2map(R,row,col)
%      map2pix()
img_utm_name='E:\2022_5_2_GF_pa1_pa2_lens_rect\pa1_012_006\mapproject\pm1_006_6087_2_12116_sub_ba_rect_utm_tif.tif';
img_utm_tfw_name='E:\2022_5_2_GF_pa1_pa2_lens_rect\pa1_012_006\mapproject\pm1_006_6087_2_12116_sub_ba_rect_utm_tif.tfw';
[X,cmap] = imread(img_utm_name);
R = worldfileread(img_utm_tfw_name,'planar',size(X));
[x,y] = pix2map(R,100,50);

1
2
3
4
5
6
7
8

(3) projfwd() 和 projinv()

实现在投影坐标系下的影像，在**像素坐标**<--->**经纬度坐标**之间的转换

%% projfwd     (forward mapprojection from lat and lon to map coordinate x,y )
img_utm_name='E:\2022_5_2_GF_pa1_pa2_lens_rect\pa1_012_006\mapproject\pm1_006_6087_2_12116_sub_ba_rect_utm_tif.tif';
img_utm_tfw_name='E:\2022_5_2_GF_pa1_pa2_lens_rect\pa1_012_006\mapproject\pm1_006_6087_2_12116_sub_ba_rect_utm_tif.tfw';

proj=geotiffinfo(img_utm_name);
[x,y]=projfwd(proj,50.2,124.3);
1
2
3
4
5
6

%% projinv    (inverse mapprojection from map coordinate to lat and lon)
img_utm_name='E:\2022_5_2_GF_pa1_pa2_lens_rect\pa1_012_006\mapproject\pm1_006_6087_2_12116_sub_ba_rect_utm_tif.tif';
img_utm_tfw_name='E:\2022_5_2_GF_pa1_pa2_lens_rect\pa1_012_006\mapproject\pm1_006_6087_2_12116_sub_ba_rect_utm_tif.tfw';

proj=geotiffinfo(img_utm_name);
[x,y]=projinv(proj,x,y);

1
2
3
4
5
6
7

(4) geotiffread()

[img,rotate_matrix]=geotiffread(filename);
[img,colormap,rotate_matrix]=geotiffread(filename);
[img,referencematrix,bbox]=geotiffread(filename);
[img,colormap,referenceematrix,bbox]=geotiffread(filename);


1
2
3
4
5
6

(5) geotiffinfo()

%% geotiffinfo   返回的info是 a structure whose fields contain  image properties and catographic infofrmation and a geotiff file

info=geotiffinfo(filename);

1
2
3
4

(6) worldfileread()

使用worldfileread这个函数需要有两个文件 **tif**和**tfw**

%% worldfileread   (read world file and return referencing object or matrix)
% 用这个函数，geotiff影像需要有2个文件： filename.tif + filename.tfw
% 主要区别是： 'geographic'--> latitude-longtitude system  
%             'planar'--> projected map coordinate
img_utm_name='E:\2022_5_2_GF_pa1_pa2_lens_rect\pa1_012_006\mapproject\pm1_006_6087_2_12116_sub_ba_rect_utm_tif.tif';
img=imread(img_utm_name);
worldfilename=getworldfilename(img_utm_name);
R=worldfileread(worldfilename,'planar',size(img));
1
2
3
4
5
6
7
8

(7) .tif 和 .tfw

%% （1）tif tfw 文件的详细介绍： https://zhuanlan.zhihu.com/p/365157364
%   （2）python为geotif文件创建tfw文件：https://qastack.cn/gis/9421/creating-tfw-and-prj-files-for-folder-of-geotiff-files
1
2

二、matlab 实现对geotiff影像的拼接

%% mosaic 2 images   (要求： 2张影像地理坐标投影：geographic，或者是  投影坐标系    ，但是两张影像行列数一致等基本信息一致)

X_west = imread('D:\matlab_mosaic\left.tif');
X_east = imread('D:\matlab_mosaic\right.tif');
X = [X_west X_east];

% 影像为地理坐标投影    下面函数的参数为： 'geographic'
% 若影像为投影坐标系， 则函数参数为： 'planar'

R_west = worldfileread('D:\matlab_mosaic\left.tfw', 'geographic', size(X_west));
R_east = worldfileread('D:\matlab_mosaic\right.tfw', 'geographic', size(X_east));
R = R_west;
R.XLimWorld = [R_west.XLimWorld(1) R_east.XLimWorld(2)];
R.RasterSize = size(X);
out_f='D:\matlab_mosaic\mosaic_li2.tif';

% 下面的coordRefSysCode 代码4326为WGS 经纬度坐标系。与原始影像对应。
% 若原始影像为投影坐标系，则修改 coordRefSysCode   即  找一些对应的EPSG code
coordRefSysCode = 4326;
geotiffwrite(out_f, X, R,'CoordRefSysCode', coordRefSysCode);


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

三. geotiff影像进行resize，并修改坐标信息

    [A,RA] = readgeoraster([files(i).folder,'\',files(i).name]);
    B = imresize(A, 0.5);

    % 修改影像的地理坐标信息
    RA.CellExtentInLatitude = RA.CellExtentInLatitude * 2;
    RA.CellExtentInLongitude = RA.CellExtentInLongitude * 2;
    RA.RasterSize = size(B);

    % 保存为geotif文件
    geotiffwrite([resize_files_folder,files(i).name], B, RA);
    

1
2
3
4
5
6
7
8
9
10
11
12
13