Filter löschen
Filter löschen

How to clip the NetCDF file (containing ASCAT soil moisture data) using a shapefile (vector)?

9 Ansichten (letzte 30 Tage)
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

Antworten (1)

Visweshwaran R
Visweshwaran R am 8 Feb. 2020
Bearbeitet: Visweshwaran R am 8 Feb. 2020
This is my data after clipping to 13-19 and 72-78 extent. But I want data exactly inside the shapefile and logical result from the 'inpolygon' is not giving properly few points are still outside. Help me with clipping lati, longi and sdata.

Community Treasure Hunt

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

Start Hunting!

Translated by