How to clip the NetCDF file (containing ASCAT soil moisture data) using a shapefile (vector)?
9 views (last 30 days)
Show older comments
The link to the shapefile is added here: Marol shapefile
The link to the NetCDF file is added here: ASCAT SSM
I tried to use inpolygon for clipping the NetCDF file but the clipping is not happening properly. I donno where my code is wrong
clear all;
close all;
clc;
%% Reading shapefile
shapefile2 = 'C:\Users\Administrator\Downloads\Programs\basin_marol\marol_wgs1984.shp';
S2 = shaperead(shapefile2);
N = length(S2); %%number of sub basins
for i = 1:N
plot(S2(i).X,S2(i).Y)
hold on
end
xs = S2(i).X(1,1:47);
ys = S2(i).Y(1,1:47);
%% Reading NetCDF file
file = 'E:\First objective\Raw data\ASCAT Data\1370844-1of1\W_XX-EUMETSAT-Darmstadt,SURFACE+SATELLITE,METOPA+ASCAT_C_EUMP_20090101040601_11429_eps_o_125_ssm_l2.nc';
% ncdisp(file)
sm = ncread(file,'soil_moisture');
latitude = ncread(file,'lat');
longitude = ncread(file,'lon');
data = cat(3,latitude,longitude,sm);
[xr yc layer] = size(sm); %%no of layers
%first layer is latitude followed by longitude and soil data in the last.
%%
k =0;
for i = 1:xr
for j = 1:yc
if(data(i,j,1)>13 && data(i,j,1)<19 && data(i,j,2)>72 && data(i,j,2)<78)
sdata1(i,j) = data(i,j,3);
% lati(i,j) = data(i,j,1);
% longi(i,j) = data(i,j,2);
lati(1,k+1) = data(i,j,1);
longi(1,k+1) = data(i,j,2);
sdata(1,k+1) = data(i,j,3);
k = k+1;
else
% sdata(i,j) = nan;
sdata1(i,j) = 0;
% lati(i,j) = nan; %%dont give 0 beacuse during plotting it will attach itself to lat = 0 and lon = 0;
% longi(i,j) =nan;
end
end
end
%% Extract data
[X Y] = meshgrid(longi,lati); %long is x and lat is y
for i =1:layer %for all the layers
for j = 1:N %for all subbasin
in = inpolygon(X,Y,S2(j).X,S2(j).Y);
end
end
%% checking weather data is correct or wrong
bb = find(in==1);
b = bb/2237; %%331 is the row of isin
col = fix(b); % is the column of isin
row = (b - col)*2237;
col = fix(b)+1;% is the row of the isin
1 Comment
KSSV
on 7 Feb 2020
the clipping is not happening properly .....why? what is the error you getting? what you want exactly?
Answers (1)
Visweshwaran R
on 8 Feb 2020
Edited: Visweshwaran R
on 8 Feb 2020
1 Comment
Sanjeev G
on 28 Aug 2020
Did you find the solution?
The problem seems to be in the subsetiing the data. Even Im in need solution.
See Also
Categories
Find more on NetCDF in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!