I have an area arbitrarily divided into contract tracts by FSA (Farm Service Agency) and an artificial waterway. I need to be able find an area to equalize the areas between two of the tracts on either side of the irregularly shaped waterway. I transferred coordinates normalized from a printed map and used FZERO to move the initial estimated boundaries of the waterway to match the FSA acreages in the tract north of the waterway, then adjusted the southern boundary to match the acreage in it. There then is a smaller second contract at the south edge that is rectangular so can easily find that boundary by ratio of the area in the tracts making up it and the total of all tracts. Having done, the remaining tract on the south side of the waterway matches in total acreage by subtraction, but doing a POLYAREA on it doesn’t match exactly, being off by roughly 0.3 A. I must be overlooking something, but I don’t see how the end result cannot match when all are calculated based on the measured total acreage.
I’ve attached the initial code and the plot of those tracts — the two little areas (tracts 5, 6) are around a well site; I didn’t try to match a few very small little bumps but matched the areas.
% terraces CRP map, return distance to match hay acreage
acre2sqft=@(A)A.*5280.^2/640; % sq ft in A acres
acre2width=@(A,L)acre2sqft(A)./(L*5280); % width in ft of A acres of tract L miles long
% Tracts=[56.23 70.65 11.38 8.65 0.33 6.56 0.43 5.24];
% isCRP=[true true true false false true false true];
% isWW=[false(1,3) true false(1,4)];
% isNorth=[true true false true false(1,4)];
% reviseed2026 map tracts/acreages
Tracts=[56.24 70.65 16.62 8.65 0.33 6.56 0.43];
isCRP=[true true true false false true false];
isWW=[false(1,3) true false(1,3)];
isNorth=[true true false true true false false];
isSouth=~isNorth;
Contract=isNorth+2*isSouth; % which contract associated with N-1, S=2
%A=sum(Tracts); % total acreage
CRP=sum(Tracts(isCRP)); % CRP acreage
WW=Tracts(isWW);
S=sum(Tracts(isSouth));
N=sum(Tracts(isNorth));
A=N+S;
toHay=sum(Tracts(isCRP & isNorth))/2; % North contract 50% target acres
% make map
ww=[0.0 3.0 4.9 7.8 8.5 8.5 7.7 5.2 3.0 0.0;
5.5 5.5 4.8 4.8 4.3 5.7 5.2 5.2 5.95 5.95]/8.5; % scaled measured dimensions
wwS=ww(:,1:5); wwN=[ww(1,1:5); fliplr(ww(2,6:end))];
% refine coordinates to match contract acreages
fn=@(dx) polyarea([wwN(1,:) 1 0], [wwN(2,:)+dx 1 1])*A-Tracts(1);
dx=fzero(fn,0);
wwN(2,:)=wwN(2,:)+dx; % adjust N boundary
% now match waterway tract area by adjustng south boundary
fn=@(dx) polyarea(ww(1,:), [wwS(2,:)+dx flip(wwN(2,:))])*A-WW;
dx=fzero(fn,0);
wwS(2,:)=wwS(2,:)+dx; % adjust S boundary
ww(2,:)=[wwS(2,:) flip(wwN(2,:))]; % update ww to match
AT1=polyarea([wwN(1,:) 1 0],[wwN(2,:) 1 1])*A;
assert(isapprox(AT1,Tracts(1)),’Tract 1 Area Mismatch’)
AWW=polyarea(ww(1,:),ww(2,:))*A; % check area
assert(isapprox(AWW,WW),’Waterway Area Mismatch’)
figure, axes, axis equal, axis tight, ylim([0 1]), xlim(ylim), box on
hold on
hA(1)=fill(ww(1,:),ww(2,:),’g’,’Tag’,’Waterway’); % the area roughly
bounds(1)=S/A; % north bound, south contract 22A +/-
hL(1)=line(xlim,[1 1]*bounds(1),’Tag’,’Divide’); % draw the line
hA(2)=fill([wwN(1,:) 1 0], [wwN(2,:) 1 1],’r’,’Tag’,’North’);
AN=polyarea([wwN(1,:) 1 0], [wwN(2,:) 1 1])*A
fill([wwS(1,:) 1 0], [wwS(2,:) bounds(1) bounds(1)],’b’)
AS=polyarea([wwS(1,:) 1 0], [wwS(2,:) bounds(1) bounds(1)])*A
AS-sum(Tracts([2 5]))
% Commented out the balance calculation — it shouldn’t move the waterway
% boundaries now — that was an earlier solution that is roughly in the
% ballpark but thought since had to update the tracts anyway would try to
% do consistent areas first — which raised the issue of not quite matching
% total area. — comment for posting — dpb
% initial values for west end of south of waterway…
%x=wwS(1,:); y=wwS(2,:); % angle segment to interpolate over
%xh=wwS(1,end-1); yh=interp1(x,y,xh); % initial use break to east
%fn=@(dx) polyarea([x(1:4) x(1,4)+dx 0], …
% [y(1:4) interp1(x,y,xh+dx,’linear’,’extrap’) interp1(x,y,xh+dx,’linear’,’extrap’)])*A + …
% Tracts(1)-toHay;
%dx=fzero(fn,0);
%yh=interp1(x,y,xh+dx); % South boundary to make 50% contract acres
%fill([x(1:end-1) xh+dx 0],[y(1:end-1) yh yh],’w’) % identify area needed south of waterway to match area
%bounds(2)=yh; % the
%AX=polyarea([x(1:end-1) xh+dx 0],[y(1:end-1) yh yh])*A; % and the area computed for check
%if abs(AX-(sum(Tracts(1:2)-Tracts(1))/2))>0.01
% error(‘Solution not found’)
%end
% add S contract and rough position of well sites
ylim([0 1]) % show whole quarter, not just large contract
WS2=Tracts(~isCRP&isSouth); % well site, contract 2 (South)
S1=Tracts(find(isCRP&isSouth,1)); % west south CRP track
xWS1=acre2width(S1,bounds(1)/2); % width of CRP tract; start of well site
dxWS1=acre2width(WS2,bounds(1)/2); % width of well site, ft
hWS1=rectangle(‘Position’,[xWS1/2640 0 dxWS1/2640 bounds(1)]);
WS1=Tracts(~(isCRP|isWW)&isNorth); % well site, contract 1 (North)
dyWS2=WS1/WS2*bounds(1); % scaled height of well site
theta=30; % rotate to NNE
dx=tand(theta)*bounds(1)*2640;
X=[xWS1 xWS1+dxWS1 xWS1+dxWS1+1.10*dx xWS1+dx]/2640;
Y=[bounds(1) bounds(1) bounds(1)+dyWS2 bounds(1)+dyWS2];
patch(X,Y,’w’)
Graph coordinates are normalized from 0-1 in both directions and so areas are fractions. What I don’t follow is why, then, the areas from polyarea() also don’t work out to add (without normalization) to unity; it seems as though if all the other areas match with it, the one last one would have to.
AS above is the blue area (Tract 2) including the diagonal well road (Tract 5).
I can make do with the discrepancy by simply using the area obtained by subtraction instead but it seems as though there ought to be a way to make things work programmatically. The 0.3 A works out as roughly a 5-ft wide strip over the half mile so it isn’t completely negligible.
acre2width(ans,0.5) % width of 1/2 mi strip to match difference
In the end, I’ll make tick marks align with the needed waypoints and convert from the fractions of the quarter mile to feet and rods (16.5 ft/rod) and finally to rod wheel counts for laying out the points on the ground. The rod wheel is 2.5 rev/rod and by fortunate circumstance, the front wheel diameter of one of the little utility tractors is an exact-enough surrogate so I can just paint a stripe across the tire and drive and count revolutions… :>)
NOTA BENE: I checked R2024b because it does use isapprox which wasn’t introduced until then, but that could be dispensed with pretty easily if wanted/needed.I have an area arbitrarily divided into contract tracts by FSA (Farm Service Agency) and an artificial waterway. I need to be able find an area to equalize the areas between two of the tracts on either side of the irregularly shaped waterway. I transferred coordinates normalized from a printed map and used FZERO to move the initial estimated boundaries of the waterway to match the FSA acreages in the tract north of the waterway, then adjusted the southern boundary to match the acreage in it. There then is a smaller second contract at the south edge that is rectangular so can easily find that boundary by ratio of the area in the tracts making up it and the total of all tracts. Having done, the remaining tract on the south side of the waterway matches in total acreage by subtraction, but doing a POLYAREA on it doesn’t match exactly, being off by roughly 0.3 A. I must be overlooking something, but I don’t see how the end result cannot match when all are calculated based on the measured total acreage.
I’ve attached the initial code and the plot of those tracts — the two little areas (tracts 5, 6) are around a well site; I didn’t try to match a few very small little bumps but matched the areas.
% terraces CRP map, return distance to match hay acreage
acre2sqft=@(A)A.*5280.^2/640; % sq ft in A acres
acre2width=@(A,L)acre2sqft(A)./(L*5280); % width in ft of A acres of tract L miles long
% Tracts=[56.23 70.65 11.38 8.65 0.33 6.56 0.43 5.24];
% isCRP=[true true true false false true false true];
% isWW=[false(1,3) true false(1,4)];
% isNorth=[true true false true false(1,4)];
% reviseed2026 map tracts/acreages
Tracts=[56.24 70.65 16.62 8.65 0.33 6.56 0.43];
isCRP=[true true true false false true false];
isWW=[false(1,3) true false(1,3)];
isNorth=[true true false true true false false];
isSouth=~isNorth;
Contract=isNorth+2*isSouth; % which contract associated with N-1, S=2
%A=sum(Tracts); % total acreage
CRP=sum(Tracts(isCRP)); % CRP acreage
WW=Tracts(isWW);
S=sum(Tracts(isSouth));
N=sum(Tracts(isNorth));
A=N+S;
toHay=sum(Tracts(isCRP & isNorth))/2; % North contract 50% target acres
% make map
ww=[0.0 3.0 4.9 7.8 8.5 8.5 7.7 5.2 3.0 0.0;
5.5 5.5 4.8 4.8 4.3 5.7 5.2 5.2 5.95 5.95]/8.5; % scaled measured dimensions
wwS=ww(:,1:5); wwN=[ww(1,1:5); fliplr(ww(2,6:end))];
% refine coordinates to match contract acreages
fn=@(dx) polyarea([wwN(1,:) 1 0], [wwN(2,:)+dx 1 1])*A-Tracts(1);
dx=fzero(fn,0);
wwN(2,:)=wwN(2,:)+dx; % adjust N boundary
% now match waterway tract area by adjustng south boundary
fn=@(dx) polyarea(ww(1,:), [wwS(2,:)+dx flip(wwN(2,:))])*A-WW;
dx=fzero(fn,0);
wwS(2,:)=wwS(2,:)+dx; % adjust S boundary
ww(2,:)=[wwS(2,:) flip(wwN(2,:))]; % update ww to match
AT1=polyarea([wwN(1,:) 1 0],[wwN(2,:) 1 1])*A;
assert(isapprox(AT1,Tracts(1)),’Tract 1 Area Mismatch’)
AWW=polyarea(ww(1,:),ww(2,:))*A; % check area
assert(isapprox(AWW,WW),’Waterway Area Mismatch’)
figure, axes, axis equal, axis tight, ylim([0 1]), xlim(ylim), box on
hold on
hA(1)=fill(ww(1,:),ww(2,:),’g’,’Tag’,’Waterway’); % the area roughly
bounds(1)=S/A; % north bound, south contract 22A +/-
hL(1)=line(xlim,[1 1]*bounds(1),’Tag’,’Divide’); % draw the line
hA(2)=fill([wwN(1,:) 1 0], [wwN(2,:) 1 1],’r’,’Tag’,’North’);
AN=polyarea([wwN(1,:) 1 0], [wwN(2,:) 1 1])*A
fill([wwS(1,:) 1 0], [wwS(2,:) bounds(1) bounds(1)],’b’)
AS=polyarea([wwS(1,:) 1 0], [wwS(2,:) bounds(1) bounds(1)])*A
AS-sum(Tracts([2 5]))
% Commented out the balance calculation — it shouldn’t move the waterway
% boundaries now — that was an earlier solution that is roughly in the
% ballpark but thought since had to update the tracts anyway would try to
% do consistent areas first — which raised the issue of not quite matching
% total area. — comment for posting — dpb
% initial values for west end of south of waterway…
%x=wwS(1,:); y=wwS(2,:); % angle segment to interpolate over
%xh=wwS(1,end-1); yh=interp1(x,y,xh); % initial use break to east
%fn=@(dx) polyarea([x(1:4) x(1,4)+dx 0], …
% [y(1:4) interp1(x,y,xh+dx,’linear’,’extrap’) interp1(x,y,xh+dx,’linear’,’extrap’)])*A + …
% Tracts(1)-toHay;
%dx=fzero(fn,0);
%yh=interp1(x,y,xh+dx); % South boundary to make 50% contract acres
%fill([x(1:end-1) xh+dx 0],[y(1:end-1) yh yh],’w’) % identify area needed south of waterway to match area
%bounds(2)=yh; % the
%AX=polyarea([x(1:end-1) xh+dx 0],[y(1:end-1) yh yh])*A; % and the area computed for check
%if abs(AX-(sum(Tracts(1:2)-Tracts(1))/2))>0.01
% error(‘Solution not found’)
%end
% add S contract and rough position of well sites
ylim([0 1]) % show whole quarter, not just large contract
WS2=Tracts(~isCRP&isSouth); % well site, contract 2 (South)
S1=Tracts(find(isCRP&isSouth,1)); % west south CRP track
xWS1=acre2width(S1,bounds(1)/2); % width of CRP tract; start of well site
dxWS1=acre2width(WS2,bounds(1)/2); % width of well site, ft
hWS1=rectangle(‘Position’,[xWS1/2640 0 dxWS1/2640 bounds(1)]);
WS1=Tracts(~(isCRP|isWW)&isNorth); % well site, contract 1 (North)
dyWS2=WS1/WS2*bounds(1); % scaled height of well site
theta=30; % rotate to NNE
dx=tand(theta)*bounds(1)*2640;
X=[xWS1 xWS1+dxWS1 xWS1+dxWS1+1.10*dx xWS1+dx]/2640;
Y=[bounds(1) bounds(1) bounds(1)+dyWS2 bounds(1)+dyWS2];
patch(X,Y,’w’)
Graph coordinates are normalized from 0-1 in both directions and so areas are fractions. What I don’t follow is why, then, the areas from polyarea() also don’t work out to add (without normalization) to unity; it seems as though if all the other areas match with it, the one last one would have to.
AS above is the blue area (Tract 2) including the diagonal well road (Tract 5).
I can make do with the discrepancy by simply using the area obtained by subtraction instead but it seems as though there ought to be a way to make things work programmatically. The 0.3 A works out as roughly a 5-ft wide strip over the half mile so it isn’t completely negligible.
acre2width(ans,0.5) % width of 1/2 mi strip to match difference
In the end, I’ll make tick marks align with the needed waypoints and convert from the fractions of the quarter mile to feet and rods (16.5 ft/rod) and finally to rod wheel counts for laying out the points on the ground. The rod wheel is 2.5 rev/rod and by fortunate circumstance, the front wheel diameter of one of the little utility tractors is an exact-enough surrogate so I can just paint a stripe across the tire and drive and count revolutions… :>)
NOTA BENE: I checked R2024b because it does use isapprox which wasn’t introduced until then, but that could be dispensed with pretty easily if wanted/needed. I have an area arbitrarily divided into contract tracts by FSA (Farm Service Agency) and an artificial waterway. I need to be able find an area to equalize the areas between two of the tracts on either side of the irregularly shaped waterway. I transferred coordinates normalized from a printed map and used FZERO to move the initial estimated boundaries of the waterway to match the FSA acreages in the tract north of the waterway, then adjusted the southern boundary to match the acreage in it. There then is a smaller second contract at the south edge that is rectangular so can easily find that boundary by ratio of the area in the tracts making up it and the total of all tracts. Having done, the remaining tract on the south side of the waterway matches in total acreage by subtraction, but doing a POLYAREA on it doesn’t match exactly, being off by roughly 0.3 A. I must be overlooking something, but I don’t see how the end result cannot match when all are calculated based on the measured total acreage.
I’ve attached the initial code and the plot of those tracts — the two little areas (tracts 5, 6) are around a well site; I didn’t try to match a few very small little bumps but matched the areas.
% terraces CRP map, return distance to match hay acreage
acre2sqft=@(A)A.*5280.^2/640; % sq ft in A acres
acre2width=@(A,L)acre2sqft(A)./(L*5280); % width in ft of A acres of tract L miles long
% Tracts=[56.23 70.65 11.38 8.65 0.33 6.56 0.43 5.24];
% isCRP=[true true true false false true false true];
% isWW=[false(1,3) true false(1,4)];
% isNorth=[true true false true false(1,4)];
% reviseed2026 map tracts/acreages
Tracts=[56.24 70.65 16.62 8.65 0.33 6.56 0.43];
isCRP=[true true true false false true false];
isWW=[false(1,3) true false(1,3)];
isNorth=[true true false true true false false];
isSouth=~isNorth;
Contract=isNorth+2*isSouth; % which contract associated with N-1, S=2
%A=sum(Tracts); % total acreage
CRP=sum(Tracts(isCRP)); % CRP acreage
WW=Tracts(isWW);
S=sum(Tracts(isSouth));
N=sum(Tracts(isNorth));
A=N+S;
toHay=sum(Tracts(isCRP & isNorth))/2; % North contract 50% target acres
% make map
ww=[0.0 3.0 4.9 7.8 8.5 8.5 7.7 5.2 3.0 0.0;
5.5 5.5 4.8 4.8 4.3 5.7 5.2 5.2 5.95 5.95]/8.5; % scaled measured dimensions
wwS=ww(:,1:5); wwN=[ww(1,1:5); fliplr(ww(2,6:end))];
% refine coordinates to match contract acreages
fn=@(dx) polyarea([wwN(1,:) 1 0], [wwN(2,:)+dx 1 1])*A-Tracts(1);
dx=fzero(fn,0);
wwN(2,:)=wwN(2,:)+dx; % adjust N boundary
% now match waterway tract area by adjustng south boundary
fn=@(dx) polyarea(ww(1,:), [wwS(2,:)+dx flip(wwN(2,:))])*A-WW;
dx=fzero(fn,0);
wwS(2,:)=wwS(2,:)+dx; % adjust S boundary
ww(2,:)=[wwS(2,:) flip(wwN(2,:))]; % update ww to match
AT1=polyarea([wwN(1,:) 1 0],[wwN(2,:) 1 1])*A;
assert(isapprox(AT1,Tracts(1)),’Tract 1 Area Mismatch’)
AWW=polyarea(ww(1,:),ww(2,:))*A; % check area
assert(isapprox(AWW,WW),’Waterway Area Mismatch’)
figure, axes, axis equal, axis tight, ylim([0 1]), xlim(ylim), box on
hold on
hA(1)=fill(ww(1,:),ww(2,:),’g’,’Tag’,’Waterway’); % the area roughly
bounds(1)=S/A; % north bound, south contract 22A +/-
hL(1)=line(xlim,[1 1]*bounds(1),’Tag’,’Divide’); % draw the line
hA(2)=fill([wwN(1,:) 1 0], [wwN(2,:) 1 1],’r’,’Tag’,’North’);
AN=polyarea([wwN(1,:) 1 0], [wwN(2,:) 1 1])*A
fill([wwS(1,:) 1 0], [wwS(2,:) bounds(1) bounds(1)],’b’)
AS=polyarea([wwS(1,:) 1 0], [wwS(2,:) bounds(1) bounds(1)])*A
AS-sum(Tracts([2 5]))
% Commented out the balance calculation — it shouldn’t move the waterway
% boundaries now — that was an earlier solution that is roughly in the
% ballpark but thought since had to update the tracts anyway would try to
% do consistent areas first — which raised the issue of not quite matching
% total area. — comment for posting — dpb
% initial values for west end of south of waterway…
%x=wwS(1,:); y=wwS(2,:); % angle segment to interpolate over
%xh=wwS(1,end-1); yh=interp1(x,y,xh); % initial use break to east
%fn=@(dx) polyarea([x(1:4) x(1,4)+dx 0], …
% [y(1:4) interp1(x,y,xh+dx,’linear’,’extrap’) interp1(x,y,xh+dx,’linear’,’extrap’)])*A + …
% Tracts(1)-toHay;
%dx=fzero(fn,0);
%yh=interp1(x,y,xh+dx); % South boundary to make 50% contract acres
%fill([x(1:end-1) xh+dx 0],[y(1:end-1) yh yh],’w’) % identify area needed south of waterway to match area
%bounds(2)=yh; % the
%AX=polyarea([x(1:end-1) xh+dx 0],[y(1:end-1) yh yh])*A; % and the area computed for check
%if abs(AX-(sum(Tracts(1:2)-Tracts(1))/2))>0.01
% error(‘Solution not found’)
%end
% add S contract and rough position of well sites
ylim([0 1]) % show whole quarter, not just large contract
WS2=Tracts(~isCRP&isSouth); % well site, contract 2 (South)
S1=Tracts(find(isCRP&isSouth,1)); % west south CRP track
xWS1=acre2width(S1,bounds(1)/2); % width of CRP tract; start of well site
dxWS1=acre2width(WS2,bounds(1)/2); % width of well site, ft
hWS1=rectangle(‘Position’,[xWS1/2640 0 dxWS1/2640 bounds(1)]);
WS1=Tracts(~(isCRP|isWW)&isNorth); % well site, contract 1 (North)
dyWS2=WS1/WS2*bounds(1); % scaled height of well site
theta=30; % rotate to NNE
dx=tand(theta)*bounds(1)*2640;
X=[xWS1 xWS1+dxWS1 xWS1+dxWS1+1.10*dx xWS1+dx]/2640;
Y=[bounds(1) bounds(1) bounds(1)+dyWS2 bounds(1)+dyWS2];
patch(X,Y,’w’)
Graph coordinates are normalized from 0-1 in both directions and so areas are fractions. What I don’t follow is why, then, the areas from polyarea() also don’t work out to add (without normalization) to unity; it seems as though if all the other areas match with it, the one last one would have to.
AS above is the blue area (Tract 2) including the diagonal well road (Tract 5).
I can make do with the discrepancy by simply using the area obtained by subtraction instead but it seems as though there ought to be a way to make things work programmatically. The 0.3 A works out as roughly a 5-ft wide strip over the half mile so it isn’t completely negligible.
acre2width(ans,0.5) % width of 1/2 mi strip to match difference
In the end, I’ll make tick marks align with the needed waypoints and convert from the fractions of the quarter mile to feet and rods (16.5 ft/rod) and finally to rod wheel counts for laying out the points on the ground. The rod wheel is 2.5 rev/rod and by fortunate circumstance, the front wheel diameter of one of the little utility tractors is an exact-enough surrogate so I can just paint a stripe across the tire and drive and count revolutions… :>)
NOTA BENE: I checked R2024b because it does use isapprox which wasn’t introduced until then, but that could be dispensed with pretty easily if wanted/needed. polyarea, fzero, sub-area normalization MATLAB Answers — New Questions