MATLAB Answers

Extracting smaller matrix from larger one, using a coordinate systrem

4 views (last 30 days)
Cian Vaugh
Cian Vaugh on 16 Apr 2021
Answered: Abhishek Gupta on 19 Apr 2021
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: -

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!