如何将地理图像中目标范围值的经纬度找出来； 如何采用基于shp文件对栅格进行分区统计?

---

前言

如何将地理坐标图像中指定的范围值的经纬度找出来？就是根据值来找位置（经纬度），我们只有tif影像，如何来找某个值范围的经纬度？如果我们有多幅影像？
如何采用基于shp文件对栅格进行分区统计，shp中有很多子区域，根据子区域进行分区统计大家都做过，GIS可以完成，但是GIS有个问题因为范围或者投影的问题会很容易报错，因此上代码强制来解决。

---

一、需要环境

matlab的地理工具箱 默认已经安装

二、代码

1.提取地理坐标图像中指定的值的经纬度

比如提取一副影像中大于1000的值的坐标经纬度：

close all
clear 
clc

tif_file = 'E:\LHASADATA\Figurecodes\data\Fig.4\370\DWELexposure.tif';
info = geotiffinfo(tif_file);
% 获取地理坐标信息
[Historicalda, R] = geotiffread(tif_file);
Ls=size(Historicalda);
lat1 = [info.CornerCoords.Lat(1) info.CornerCoords.Lat(2); info.CornerCoords.Lat(4) info.CornerCoords.Lat(3)];
lon1 = [info.CornerCoords.Lon(1) info.CornerCoords.Lon(2); info.CornerCoords.Lon(4) info.CornerCoords.Lon(3)];
[Xi,Yi] = meshgrid(1:2,1:2);
[XI,YI] = meshgrid(1:1/(Ls(2)-1):2,1:1/(Ls(1)-1):2);
LatLsatf = interp2(Xi,Yi,lat1,XI,YI);
LonLsatf = interp2(Xi,Yi,lon1,XI,YI);
% 获取大于1000的像素
positive_values = Historicalda> 1000;
eventcount=Historicalda(positive_values);
lat = LatLsatf(positive_values);
lon = LonLsatf(positive_values);
shpdata=[lat,lon,eventcount];

2.对栅格进行分区统计

比如要对一个shape文件（不规则）的全球国家区域的粮食产量，粮食产量是栅格数据，而且还要将最终结果导出来

clear
clc
shpfile = shaperead('G:\2024codes\folium_documentation-main\世界国家.shp');
msweptif_file='2001-2020.tif';
ssp370tif_file='2021-2100.tif';
ssp585tif_file='2021-2100.tif';
% lon 为原数据的lon
% lat 为原数据的lat
% shp_file为要裁剪的shp file
% lonnew=min(shpfile(kk).X):0.1:max(shpfile(kk).X);
% latnew=min(shpfile(kk).Y):0.1:max(shpfile(kk).Y);
% datanew=nan(numel(latnew),numel(lonnew));

mswepexposure=[];
ssp370exposure=[];
ssp585exposure=[];
for kk = 1:length(shpfile)

out = [shpfile(kk).X; shpfile(kk).Y]';
[mswepdata, R1] = geotiffread(msweptif_file);
[ssp370data, R2] = geotiffread(ssp1tif_file);
[ssp585data, R3] = geotiffread(ssp2tif_file);

lat=min(R1.LatitudeLimits):0.1:max(R1.LatitudeLimits);
lon=min(R1.LongitudeLimits):0.1:max(R1.LongitudeLimits);
[xx, yy] = meshgrid(lon, lat);
xx(:,end)=[];
xx(end,:)=[];
yy(:,end)=[];
yy(end,:)=[];
index = inpolygon(xx, yy, out(:,1), out(:,2));
index=flipud(index);%有的不需要翻转，需要提前取一个来判断
    % 将非多边形内的值设置为NaN
%data(~index) = nan; 
mswepdata(index == 0) = nan; % 将非多边形内的值设置为NaN
ssp1ata(index == 0) = nan; 
ssp2data(index == 0) = nan; 
%subdwel=data(index);
mswepexposurecount=nansum(nansum(mswepdata));
ssp370exposurecount=nansum(nansum(ssp1data));
ssp585exposurecount=nansum(nansum(ssp2data));
mswepexposure=[mswepexposure;mswepexposurecount];
ssp370exposure=[ssp370exposure;ssp1exposurecount];
ssp585exposure=[ssp585exposure;ssp2exposurecount];
%Rnew= georasterref('RasterSize', size(),'Latlim', [double(min(lat)) double(max(lat))], 'Lonlim', [double(min(lon)) double(max(lon))]);
%geotiffwrite(fullfile(['\test.tif']),data,Rnew);
disp(kk)
end
disp('done')

---

总结

记录科研中处理数据通用的代码，顺手发到CSDN上造福广大科研民工，欢迎关注我的公众号，获取更多科研代码和前沿论文资讯等相关内容