matlab与geotiff影像的函数,用法介绍
文章目录前言一、geotiffread,geotiffwrite,geotiffinfo二、matlab 实现对影像的拼接总结前言matlab给影像加地理坐标信息,一直没有很系统的明白,现在学习并记录吧一、geotiffread,geotiffwrite,geotiffinfofclose allclear allclcfilename='D:\sift_test\landsat_5.tif';o
·
文章目录
前言
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
(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);
(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);
%% 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);
(4) geotiffread()
[img,rotate_matrix]=geotiffread(filename);
[img,colormap,rotate_matrix]=geotiffread(filename);
[img,referencematrix,bbox]=geotiffread(filename);
[img,colormap,referenceematrix,bbox]=geotiffread(filename);
(5) geotiffinfo()
%% geotiffinfo 返回的info是 a structure whose fields contain image properties and catographic infofrmation and a geotiff file
info=geotiffinfo(filename);
(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));
(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
二、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);
三. 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);
更多推荐
所有评论(0)