Extracting smaller matrix from larger one, using a coordinate systrem

7 views (last 30 days)
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

Answers (1)

Abhishek Gupta
Abhishek Gupta on 19 Apr 2021
Hi,
As per my understanding, you want to extract a part of a matrix. Referring to the following MATLAB Answer, which might help you resolve the issue: -
For more information, check out the following documentation link: -

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!