what's wrong with this code? polygon function

BN
2019 11 4
1 个回答
回答已采纳
更新时间：2019 11 5
4 次查看（30 天）

0 个投票

I have a netcdf file containg data for a region over 4 country. I only want to have the data over specific country and not the rest of the region. I read so many document but still I cant do that and need help. here what's I've done yet:
filename= 'chirps-v2.0.monthly.nc'
time= ncread(filename,'time') %read time
lat = ncread(filename,'latitude'); %reading latitude
lon = ncread(filename,'longitude'); %reading longitude
precip = ncread(filename,'precip'); %reading the main variable precip=(lon*lat*time)
precip_mean= mean(precip, 3) %average of precip in all times
[x, y] = borders('Iran Islamic Republic of');
in = inpolygon(x,y,lon,lat) %!!!!! I GOT ERROR HERE !!!!! I write rest of code based on examples in MathWorks
precip_mean=precip_mean';
precip_mean(~in)=NaN;
surf(lon, lat, Mean_precip(:,:,:).'); view(2)
shading interp
colormap jet
hold on
pgon = polyshape (x,y,'simplify',false);
plot(pgon,'FaceColor','none','edgecolor','k','facealpha',1);
I want only precip that are in my study region's border.

4 个评论
显示 2更早的评论 隐藏 2更早的评论

BN
BN 2019-11-4
编辑：BN 2019-11-4
Oh sorry, I forgot that:
Matrix dimensions must agree.
Error in inpolygon>close_loops (line 233)
if ~any(xnan | ynan)
Error in inpolygon (line 83)
[xv, yv] = close_loops(xv, yv);
Walter Roberson
Walter Roberson 2019-11-4
What shows up for size() of x and y?
BN
BN 2019-11-4
编辑：BN 2019-11-4
size(x)
ans =
1 2201
size (y)
ans =
1 2201
size(lat)
ans =
360 1
size(lon)
ans =
720 1
size (precip)
ans =
720 360 1512
size(precip_mean)
ans =
720 360

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Walter Roberson
Walter Roberson 2019-11-4

2 个投票

inpolygon requires the first two arguments (query locations) to be the same size as each other, and the polygon definition in the third and fourth argument to be vectors the same size as each other. You cannot define the boundaries with 360 lat and 720 long.

7 个评论
显示 5更早的评论 隐藏 5更早的评论

BN
BN 2019-11-4
Ok thank you for your answer, I try limit lat and lon arround study region.
filename= 'chirps-v2.0.monthly.nc'
time= ncread(filename,'time') %read time
lat = ncread(filename,'latitude'); %reading latitude
lon = ncread(filename,'longitude'); %reading longitude
precip = ncread(filename,'precip'); %reading the main variable precip=(lon*lat*time)
precip_mean= mean(precip, 3) %average of precip in all times
[x, y] = borders('Iran Islamic Republic of');
ind1 = find(lon>=44 & lon<=64); %clip country
ind2 = find(lat>=24 & lat<=40);
lat=lat(ind2); %new lat and lon
lon=lon(ind1);
Now I have size(lat)=32*1 and size(lon)= 40*1
size of x and y are the same and 1*2201
I think I should do something on lat and lon like combine them in order to have both of them in the same sizes. then:
in = inpolygon(lon,lat,x,y)
how to do this?
Thank you
Walter Roberson
Walter Roberson 2019-11-4
filename= 'chirps-v2.0.monthly.nc'
time= ncread(filename,'time') %read time
lat = ncread(filename,'latitude'); %reading latitude
lon = ncread(filename,'longitude'); %reading longitude
precip = ncread(filename,'precip'); %reading the main variable precip=(lon*lat*time)
precip_mean= mean(precip, 3) %average of precip in all times
[x, y] = borders('Iran Islamic Republic of');
[lonG, latG] = ndgrid(lon, lat); %you might have to reverse these
in = inpolygon(lonG, latG, x, y);
inmask = repmat(in, 1, 1, size(precip,3));
masked_precip = precip;
masked_precip(~inmask) = nan;
mean_masked_precip = mean(masked_precip, 3, 'omitnan');
BN
BN 2019-11-5
I really appreciate your kindness and time, just one thing exists that I don't understand. after:
surf(lon, lat, mean_masked_precip(:,:,:).'); view(2)
axis xy
% [Loni2,Lati2] = meshgrid(lon,lat);
%shading interp
cmocean 'rain' % rainy to dry colormap
xlabel longitude
ylabel latitude
% [Lati2,Loni2,newpreceip,Mean_precip] = recenter(Lati2,Loni2,newpreceip,Mean_precip);
hold on
borders('countries','color',rgb('dark gray'))
cmocean 'rain'
cb = colorbar;
cb.Label.String = 'Average Precipitation (mm/m)';
xlabel longitude
ylabel latitude
title ('Precipitation (Monthly Average from 1982 to 2015)')
the output is :
Capture3.JPG
After drawing it should have fallen to the yellow highlight but it has fallen somewhere else. What's the problem?
Thank you again
Walter Roberson
Walter Roberson 2019-11-5
编辑：Walter Roberson 2019-11-5
That looks to me to be classic reversal of x and y, but I cannot see at the moment where that is happening. lon = x, lat = y
BN
BN 2019-11-5
I get it from here: https://www.chadagreene.com/CDT/CDT_Getting_Started.html
BN
BN 2019-11-5
I was trying some edits this day from morning till now but surrender. I don't know what to do. I tried it on 4 different NetCDF from different climate models. your help is really precious for me. thank you for your time.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Geographic Plots 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by