Extracting smaller matrix from larger one, using a coordinate systrem

5 vues (au cours des 30 derniers jours)
Cian Vaugh
Cian Vaugh le 16 Avr 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

Réponses (1)

Abhishek Gupta
Abhishek Gupta le 19 Avr 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: -

Catégories

En savoir plus sur Creating and Concatenating Matrices dans Help Center et File Exchange

Produits


Version

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by