Extracting smaller matrix from larger one, using a coordinate systrem

Question

Cian Vaugh 2021-4-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/803836-extracting-smaller-matrix-from-larger-one-using-a-coordinate-systrem

回答： Abhishek Gupta 2021-4-19

Hi,

Im dealing with NIMROD weater data for a project . i have adapted and modifyed the program to extract the data into a large rain intensity matrix( rr_dat_mat), i only seed a subsection of the data between a certin coo ordinates but im having isssures going between easting northings and the index's

below is the the code that im running i cant provide the data files, However hopefull this is enbough to help with with my blunders

%% GLobal Var
utmzone = '29 U' %Universal Transverse Mercator coordinate system irelnad refrance 
%%
MASTERFILE = ...
'metoffice-c-band-rain-radar_uk_20210129_1km-composite.dat.gz.tar'
%'metoffice-c-band-rain-radar_europe_20210211_5km-composite.dat.gz.tar'
%%file handeling __________________________________________________________
scrsz = get(0,'ScreenSize');
fh1 = figure('OuterPosition',[1 scrsz(4)*.05 scrsz(3)*0.5 scrsz(4)* 0.95]);
disp('Extracting files')
files = gunzip(untar(MASTERFILE,'TempGz'),'TempDat');
%ectractedGz  = gunzip('metoffice-c-band-rain-radar_uk_20210129_1km-composite.dat.gz')
NoFiles = length(files)
disp('files extracted')
%% Easting northing Isioation area _________________________________________ 
iosilationAreapoint1.easting = 156534;
iosilationAreapoint1.northing =187799 ;
AREA_LENGTH = 20;
AREA_WIDTH = 10; 
[isolationAreaPoint1.Lat,isolationAreaPoint1.Lon] = ...
utm2deg(iosilationAreapoint1.easting,iosilationAreapoint1.northing,utmzone);
iosilationAreapoint2.easting =171582;
iosilationAreapoint2.northing = 166267;
[isolationAreaPoint2.Lat,isolationAreaPoint2.Lon] = ...
utm2deg(iosilationAreapoint2.easting,iosilationAreapoint2.northing,utmzone);
%geoplot([isolationAreaPoint2.Lat isolationAreaPoint1.Lat],...
 %   [isolationAreaPoint2.Lon isolationAreaPoint1.Lon],'g-*')
hold on 
%% Frame loop itterator
for i = 1:6:NoFiles
% run rdnim1k to extract data items form .dat file 
[int_gen_hd, rl_gen_hd, rl_datsp_hd, char_hd, int_datsp_hd, ...
rr_dat_mat] = rdnim1km( files{i} );
%test translation factor 
Trans.x = rl_datsp_hd(2)-1; %matlab indexes start at 1 
Trans.y = rl_datsp_hd(1)- 1;
%Get co ordnates for bottom left
grahOrigin.x = rl_datsp_hd(2);
grahOrigin.y = rl_datsp_hd(5);
%used to identify top right corener ( make sure orintation and endings
%configured correctly 
%     xmax = int_gen_hd(19);
%     xmax = 157363.75;
%     rr_dat_mat(1:20,xmax-19:xmax) = 600;
% Plot rain radar image on National Grid (NG) scale
% -------------------------------------------------
% Set colour limits for image and colourbar (1000=31.25mm/hr - you may
% need to increase this for extreme weather events).
clims = [ 0 1000 ];
% create x and y vectors one value for each pixel along edge of plot
%NB: yg is deliberately in descending order to allow correct display of
% both array and NG co-ordinates on plot y-axis
xg = linspace(rl_datsp_hd(2),rl_datsp_hd(4),int_gen_hd(19));
yg = linspace(rl_datsp_hd(1),rl_datsp_hd(5),int_gen_hd(18));
% display the rain radar data image
%test refrence system
rr_dat_mat(rl_datsp_hd(2)- Trans.x,rl_datsp_hd(1)-Trans.y)= 600;
%index and save rowes fo selected area 
M = zeros;
for eastSave = (iosilationAreapoint1.easting- Trans.x):(iosilationAreapoint2.easting-Trans.x)
    for northSave = (iosilationAreapoint1.northing- Trans.y):( iosilationAreapoint2.northing-Trans.y)
        M = [M  rr_dat_mat(eastSave,northSave)];
        %highlight Area on Whole Pitcher 
        rr_dat_mat(eastSave,northSave) = 600; 
    end
    %dlmwrite('test.csv',M,'delimiter',',','-append');
end
%plot data
imagesc(xg,yg,rr_dat_mat,clims);
% use rainbow colour scheme
colormap(jet);
% display axis scales corresponding to...
% NG co-ordinates of area covered by image: left right bottom top 
wholePicture = [rl_datsp_hd(2) rl_datsp_hd(4) rl_datsp_hd(5) rl_datsp_hd(1) ];%aquiered form header file 
%windowedPitcher = [157363.75 172053.28 165807.68 187324.25]; % easting norting of bottom left and top right 
windowedPitcher = [iosilationAreapoint1.easting iosilationAreapoint2.easting iosilationAreapoint2.northing iosilationAreapoint1.northing];
%comment or oncomment one of the 2 following lines in order to get differnt
%views 
axis(wholePicture);
%axis(windowedPitcher);
axis xy
% ensure square grid i.e 1km easting is same as 1km northing
axis equal
% display uncompressed filename (which includes timestamp) as title
title('frame');
% display key to colours on rhs of plot
colorbar;
% allow successive plots to use same figure and plot area
hold on
grid on 
% mark location of interest between white and red 'x'
plot(iosilationAreapoint1.easting,iosilationAreapoint1.northing,'xw');
hold on 
plot(iosilationAreapoint2.easting,iosilationAreapoint2.northing,'xr');
disp( 'plotting seriies - loop number :' );
disp(i)
drawnow
pause(0.1)
end