circle problem, circumscribed to another circle and the area, area per image detection

I have been working on a transportation problem, I have a series of coordinates that always vary, but for this question let's suppose that I have 13 coordinates, each chord represents the center of a circle with a radius of 400m.

the circles are on an x, y plane. The idea is to find the total area that is inside each circle, the problem is that it is not simply multiplying the number of circles by (PI * R ^ 2), because many of them have 1 or more shared areas that do not they can be added more than once.

These circles may or may not intercept each other, through the program that I am designing, I even obtain the distance matrix between each one, at first I thought that the best way to solve the problem was mathematically, but today I realized that it is a problem of mathematical complexity .

So looking on the internet I found that more people have had a similar problem and that it can be solved analytically.

It turns out that there are people who say that the problem could be solved if each point is plotted on a 2d plane, then a 400 meter Buffet is placed at each point, each circle is painted in one color and then the background in another color.

Finally, an algorithm is required that is capable of finding the area that is colored, suppose that what is inside each circle is yellow and the outside is green.

could you please help me with this

attached document excel file

1)Once I have plotted the figure, how do I put the 400 meter buffer at each point, I know the x, y coordinate of each point)

2)Supposing that I manage to obtain a graph with the points and their respective buffer of 400 meters of radius, how I paint everything that is inside the buffer of the same color

3)how do i find the colored area inside each buffer

T=readtable('mini2.xlsx')
H = table(T.lon,T.lat);
plot(T.lon,T.lat,'ko*400') 

plot It has no buffer, I don't know how to implement it

Thanks for your time.

examples:

I have read about the problem here : https://qastack.mx/programming/1667310/combined-area-of-overlapping-circles

http://jsfiddle.net/broofa/WRpmW/5/

[KSSV]

The give centers are in degrees right? You have to convert them into meteres first right?

[David Alejandro Ramirez Cajigas]

the centers if I have them, in a coordinate file

[Andrew]

Are approximate answers OK? Since the radii are always 400 m, I would consider finding out ahead of time how much two circles overlap as a function of distance. At x/r~=0, the second circle has basically zero area, and at x/r~=2, it's basically a whole circle. Then you can simply pairwise see which ones overlap and apply the function. The only serious issue I see is that if you have three or more overlapping, you'll overcount the overlap, so I suppose this would work best if it's a pretty sparse map or the little bit in the midst of three circles doesn't matter much.

Answer 1

[DGM]

This could be improved or elaborated, but it's a simple start:

% let's assume that 1px = 1 square meter
r = 400;
C = rand(10,2)*5000 + 1000; % test locations (x,y)

% get padded image size
p = 10; % extra padding
s = ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
C = C - [min(C(:,1)) min(C(:,2))] + r + p; % shift image

% draw all circles
m = false(s);
for c = 1:size(C,1)
	thisroi = images.roi.Circle('center',C(c,:),'radius',r);
	m = m | createMask(thisroi,m);
end

totalarea = sum(m(:)) % total number of white pixels = area

totalarea = 4914749

imshow(m)

There is also bwarea(), but you might want to read the docs to make sure it's what you want. It's not the same as a simple sum.

Of course, if there are many such circles and they are very spread out, the image may be very large. It might be rather slow to do it by drawing circles like this. This is a different approach for the same thing, which should be faster than explicit circle drawing.

% same stuff
r = 400;
C = rand(10,2)*5000 + 1000; % test locations (x,y)

% same stuff
p = 10; 
s = ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
C = round(C - [min(C(:,1)) min(C(:,2))] + r + p);

% use distance xform
m = false(s);
m(sub2ind(s,C(:,2),C(:,1))) = 1;
m = bwdist(m)<=r;

totalarea = sum(m(:)) % total number of white pixels = area

totalarea = 4164733

imshow(m)

[David Alejandro Ramirez Cajigas]

Can I put that only the white shaded area is being added?

the program gives me an error.

Thanks a lot

T=readtable('mini2.xlsx')

C = table(T.lon,T.lat);
r=400
p=10
s=ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
C = round(C - [min(C(:,1)) min(C(:,2))] + r + p);
m = false(s);
m(sub2ind(s,C(:,2),C(:,1))) = 1;
m = bwdist(m)<=r;
totalarea = sum(m(:))
imshow(m)

[David Alejandro Ramirez Cajigas]

Thank you

the excel has real coordinates of a random area of the planet.

your command are random values, so both are equal.

but I do not know the reason why when I replace by the call of the excel list it does not take it as the same

This command throws me an Array, table or a Cell? ok now i know the command is only .

nametable{:,:}

T=readtable('mini2.xlsx')
C = table(T.lon,T.lat);
C = C{:,:}

but.... i have this error

C = rand(10,2)*5000 + 1000

   1.0e+03 *

    4.2787    4.5302
    1.1786    1.1592
    5.2456    2.3846
    5.6700    1.2309
    4.3937    1.4857
    4.7887    5.1173
    4.7157    4.4741
    2.9611    2.5855
    4.2774    5.7511
    1.8559    1.1722

The one that I use brings a table from an excel and apparently it is not compatible with the composition.

T=readtable('mini2.xlsx')
C = table(T.lon,T.lat);

13×2 table

     Var1       Var2 
    _______    ______

    -76.498    3.4309
    -76.502    3.4321
    -76.501    3.4356
      -76.5    3.4355
    -76.502     3.432
    -76.505    3.4336
    -76.506    3.4336
    -76.503     3.434
    -76.501    3.4328
    -76.503    3.4337
    -76.501    3.4328
    -76.499    3.4318
    -76.498    3.4314

[DGM]

Your data is a table (instead of a numeric array) and the values are in degrees (instead of meters). You would need to convert the table to a numeric array

numericarray = table2array(mytable);

To convert from lat/lon to meters, you may be able to use one of the following

https://www.mathworks.com/help/aerotbx/ug/lla2flat.html

https://www.mathworks.com/help/driving/ref/latlon2local.html

https://www.mathworks.com/matlabcentral/answers/44064-how-would-i-convert-latitude-and-longitude-to-meters-from-a-point

but you may have to pay attention to how the results are shifted, if a reference point is used.

You will also need to keep track of which column is x and which is y.

[David Alejandro Ramirez Cajigas]

the funtion latlon2local it does not work in any way, it helps, I have even entered the example from the documentation. even so it still does not work.

latlon2local(42.3648,-71.0214,0) %Neither the example in the documentation works, I am using matlab 2021a

the axes convention

x= longitud

y= latitud

How can I find the relationship between pixel and metric unit of measure.

for example if I have a 1 meter by 1 meter square which is a meter squared, how many pixels are there?

so I could have the total metric relationship based on your algorithm. Thank you

T=readtable('mini2.xlsx')
b = table(T.lon,T.lat);
C = b{:,:};
r=400
p=10
s=ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
C = round(C - [min(C(:,1)) min(C(:,2))] + r + p);
m = false(s);
m(sub2ind(s,C(:,2),C(:,1))) = 1;
m = bwdist(m)<=r;
totalarea = sum(m(:))
imshow(m)

[DGM]

I don't have the toolbox for latlon2local(). I don't know if that's the problem you're having. I used lla2flat(), but that's part of a different toolbox which you might also not have.

% recreate example table
lon = [-76.498 -76.502 -76.501 -76.5 -76.502 -76.505 -76.506 -76.503 -76.501 -76.503 -76.501 -76.499 -76.498];
lat = [3.4309 3.4321 3.4356 3.4355 3.432 3.4336 3.4336 3.434 3.4328 3.4337 3.4328 3.4318 3.4314];
lonlat = table(lon.',lat.');

% process table and convert to meters
lla = [fliplr(table2array(lonlat)) zeros(size(lonlat,1),1)];
position = lla2flat(lla,[min(lla(:,1)) min(lla(:,2))],90,0);

% same as before
r = 400;
C = position(:,1:2); % locations (x,y)

p = 10; % padding
s = ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
C = round(C - [min(C(:,1)) min(C(:,2))] + r + p);

m = false(s);
m(sub2ind(s,C(:,2),C(:,1))) = 1;
m = bwdist(m)<=r;

totalarea = sum(m(:)) % total number of white pixels = area

totalarea = 1555013

imshow(m)

The points are a lot closer than I expected, but that seems about right.

If you don't have lla2flat(), then there are also user-submitted tools. I just did a quick search and found this:

https://www.mathworks.com/matlabcentral/fileexchange/45699-ll2utm-and-utm2ll

which gives similar results.

% process table and convert to meters
lonlat = fliplr(table2array(lonlat));
position = ll2utm(lonlat(:,1),lonlat(:,2)); % [xpos, ypos] in meters
% and so on...

[David Alejandro Ramirez Cajigas]

1) I am very grateful to you

2) Then could you suppose that total area = 1555013 is the total in square meters of the problem or the number of pixels?

3) r = 400; its the same like a 400 meters

3) I add the real image of the views in google

THANKS A LOT

[DGM]

Well, since the lla conversions output coordinates in meters, and the convention i used is that image coordinates are in the same units as C, then 1px = 1 meter. The consequence of that is then yes, r = 400px = 400m and the number of white pixels in the output is the same as area in square meters.

[David Alejandro Ramirez Cajigas]

Thanks

[David Alejandro Ramirez Cajigas]

i have more one question

[Image Analyst]

@David Alejandro Ramirez Cajigas, ok, fine, but you forgot to ask it. What is it?

[David Alejandro Ramirez Cajigas]

It is posted in this same forum below

[David Alejandro Ramirez Cajigas]

I have posted the question right here at the end, it will be very kind if you could help me

Thank you

down in this forum, is the question

Thank you

Answer 2

[David Alejandro Ramirez Cajigas]

for Java I found this on the web could matlab be adapted?

author : http://jsfiddle.net/broofa/WRpmW/5/

// Get HTML canvas element (and context) to draw into
var canvas = document.getElementById('canvas');
var ctx = canvas.getContext('2d');

// Lil' circle drawing utility
function circle(x,y,r) {
  ctx.beginPath();
  ctx.arc(x, y, r, 0, Math.PI*2);
  ctx.fill();
}

// Clear canvas (to black)
ctx.fillStyle = 'black';
ctx.fillRect(0, 0, canvas.width, canvas.height);

// Fill shape (in white)
ctx.fillStyle = 'white';
circle(40, 50, 40);
circle(40, 10, 10);
circle(25, 15, 12);
circle(35, 90, 10);

// Get bitmap data
var id = ctx.getImageData(0, 0, canvas.width, canvas.height);
var pixels = id.data; // Flat array of RGBA bytes

// Determine area by counting the white pixels
for (var i = 0, area = 0; i < pixels.length; i += 4) {
  area += pixels[i]; // Red channel (same as green and blue channels)
}

// Normalize by the max white value of 255
area /= 255;

// Output result
document.getElementById('result').innerHTML = area.toFixed(2);

[DGM]

This is JS, and is conceptually equivalent to the first example.

[David Alejandro Ramirez Cajigas]

Thank you for all

Answer 3

[David Alejandro Ramirez Cajigas]

@DGM @Image Analyst

Will it be possible to make the plotted image of black and white circles to measure the area, be drawn within specific margins? that is to say that it is drawn, for example, inside a straight line and where a circle touches the "wall" of the rectangle, is that area not counted?

I could have the coordinates of the rectangle and from there I would already have the "wall" where the circuit cannot advance.I made a small drawing and I will post it within the same link as the question already solved. with an exampleThank you.

1) the point cloud in longitude and latidut is

lat	lon
3.4310571	-76.5249196
3.4330517	-76.5281573
3.4329578	-76.5278992
3.4259733	-76.5169571
3.4141012	-76.5322045
3.4132041	-76.5212801
3.4147463	-76.5212291
3.4200984	-76.5090573
3.4193326	-76.5304345
3.4191612	-76.5226856
3.4195618	-76.5180319
3.4194799	-76.5200061
3.4197733	-76.514478
3.4199018	-76.5121552
3.4191972	-76.5076382
3.4186178	-76.5074345
3.4172711	-76.5061631
3.4174853	-76.5059834
3.4131446	-76.5021961
3.413158	-76.5019253
3.410971	-76.5001202
3.4120227	-76.4916586
3.4146176	-76.4889863
3.4187462	-76.5249977
3.4121893	-76.5258832
3.4266255	-76.5285806
3.4291025	-76.5167687
3.4347968	-76.5141764
3.4321028	-76.5101736
3.4271533	-76.5053988
3.4185984	-76.4869869
3.4218167	-76.4909459
3.4248689	-76.4945159
3.4282086	-76.5327531
3.4354939	-76.5187478
3.4334187	-76.5101116
3.4352097	-76.5118256
3.4277076	-76.5368692
3.4309797	-76.5358927
3.4333837	-76.5351849
3.4349393	-76.5347021
3.4327974	-76.5364751
3.4182723	-76.5367191
3.4090406	-76.5355765
3.4077261	-76.5315237
3.4109503	-76.4893841
3.4086604	-76.5217632
3.4107276	-76.5215164
3.4166607	-76.5211731
3.4197919	-76.520905
3.4226876	-76.5205857
3.4247333	-76.5203818
3.4252902	-76.5180966
3.4295177	-76.5165489
3.4295177	-76.5165489
3.431641	-76.5163424
3.434254	-76.5147331
3.4347146	-76.5147384
3.4324815	-76.5163367
3.4213276	-76.517367
3.4175283	-76.5222486
3.412915	-76.5227287
3.410417	-76.523005
3.4071183	-76.5233778
3.4253297	-76.5015155
3.4339325	-76.4964618
3.4309045	-76.4980119
3.4244061	-76.5032478
3.4351374	-76.5232265
3.4346636	-76.5230613
3.4352634	-76.5206634
3.4330464	-76.5264623
3.4329822	-76.525545
3.432503	-76.5272398
3.431776	-76.528092
3.4305833	-76.5299112
3.4284118	-76.5320919
3.42642	-76.5337575
3.4241897	-76.533956
3.4235712	-76.5346507
3.4203582	-76.5356673
3.4178333	-76.5365309
3.4186983	-76.5357638
3.4215123	-76.5348599
3.4260317	-76.5333472
3.431708	-76.53415
3.4312579	-76.5281195
3.4321469	-76.5268189
3.4104115	-76.5293914
3.4129524	-76.5282836
3.4116382	-76.5231909
3.4225161	-76.5052968
3.4223848	-76.5046797
3.4201063	-76.5077996
3.4198948	-76.5072551
3.4282383	-76.4985591
3.4307088	-76.4971941
3.4295397	-76.498683
3.4341627	-76.4954614
3.4186526	-76.5283958
3.420709	-76.5345703
3.4148539	-76.5343228
3.4143741	-76.5182786
3.4143097	-76.5204408
3.4127084	-76.5180348
3.4121921	-76.5152421
3.4101545	-76.5174233
3.4144032	-76.5129681
3.4171883	-76.5100628
3.4339407	-76.53169
3.4313357	-76.5324383
3.435523	-76.5321513
3.4296035	-76.527203
3.4293248	-76.5276801
3.4261009	-76.5289884
3.4230866	-76.5299733
3.4235895	-76.5295651
3.4198576	-76.5310217
3.4198307	-76.5307561
3.4168025	-76.5317512
3.416891	-76.5319739
3.4140367	-76.5326686
3.4140957	-76.5328912
3.4119804	-76.5336127
3.4118654	-76.5333713
3.4093833	-76.5341652
3.4096271	-76.534353
3.4072306	-76.5351255
3.4085317	-76.5191425
3.4089332	-76.5194188
3.4076618	-76.5207116
3.4112626	-76.5170021
3.4133751	-76.5148108
3.4176403	-76.510361
3.4163313	-76.5116723
3.4153969	-76.5043551
3.4084892	-76.4974891
3.4101436	-76.4990257
3.4144384	-76.503151
3.4351054	-76.5284981
3.4261281	-76.5275993
3.4261523	-76.527849
3.4229822	-76.5274654
3.4228698	-76.5277042
3.4193409	-76.5272482
3.4195142	-76.5274917
3.4174159	-76.5276934
3.4175039	-76.5274625
3.4152952	-76.5282755
3.4153355	-76.5280341
3.4136218	-76.5285115
3.4128266	-76.5289863
3.4105909	-76.5295951
3.4122988	-76.4918565
3.4172202	-76.4868609
3.414781	-76.4891394
3.4114073	-76.4938627
3.4161733	-76.4946323
3.413493	-76.4958315
3.4153673	-76.4960917
3.409121	-76.5029153
3.4107622	-76.5010752
3.4078625	-76.5035831
3.4223742	-76.5071477
3.4225534	-76.5073462
3.4296915	-76.5034007
3.4300717	-76.5034999
3.427595	-76.504774
3.4321119	-76.5024404
3.4273729	-76.5045942
3.4354533	-76.5004663
3.4320209	-76.5022017
3.4252149	-76.5060265
3.4250836	-76.5057395
3.4355899	-76.5261781
3.4306315	-76.5255236
3.4283743	-76.5244554
3.4276114	-76.5261325
3.4235978	-76.5259367
3.4216031	-76.5252152
3.4217371	-76.5227609
3.4199565	-76.5219831
3.4077286	-76.5265107
3.4108104	-76.525537
3.42046	-76.5240967
3.423772	-76.5242335
3.4260048	-76.5243488
3.4275149	-76.5116888
3.4254534	-76.5118604
3.42231	-76.5121072
3.4119459	-76.5144917
3.4120689	-76.5143116
3.410232	-76.5128663
3.4103927	-76.5127
3.4080554	-76.5104764
3.4079321	-76.5106856
3.410584	-76.5180612
3.4136912	-76.5177206
3.4146069	-76.519405
3.4281973	-76.5199898
3.4324517	-76.5213684
3.4319432	-76.5237341
3.4256217	-76.5249679
3.4257289	-76.5229053
3.4234424	-76.5215964
3.4102992	-76.5182481
3.4148131	-76.5178014
3.4211585	-76.517155
3.414896	-76.5176163
3.4173218	-76.5175815
3.4173192	-76.517375
3.4198733	-76.5172382
3.4197957	-76.5174501
3.4226928	-76.517037
3.4080927	-76.5165086
3.4072227	-76.5159668
3.4076885	-76.5163987
3.4144512	-76.4930341
3.4129628	-76.4915964
3.4176748	-76.4959094
3.419292	-76.4974544
3.4224996	-76.5005174
3.4207702	-76.4988491
3.4227033	-76.5040283
3.4205988	-76.5019228
3.4180153	-76.4994632
3.4162721	-76.4980081
3.4352152	-76.5141913
3.4335874	-76.5092426
3.4330037	-76.5074616
3.433582	-76.5052032
3.4336247	-76.5059327
3.4339756	-76.5031808
3.4327895	-76.501204
3.4337159	-76.5030869
3.4328028	-76.5008097
3.4318442	-76.4993533
3.4314239	-76.4980551
3.4344601	-76.5118336
3.4192084	-76.5207168
3.4193369	-76.5181929
3.4194627	-76.5159237
3.4196073	-76.5133407
3.4197599	-76.5103903
3.428547	-76.5330845
3.4121579	-76.4873705
3.4297281	-76.5109285
3.4261112	-76.5168337
3.4229298	-76.5171669

2) the result image is

3) the idea is that the points are plotted in a "rectangular" range demarcated by coordinates

for example:

lat         lon
3.4356   -76.5369  %lan and lot 1
3.4071    -76.4868 %lan and lot 2

How do I get it to plot and measure me within the geometric limits of the rectangle, formed by the given coordinates?

Thank you

[David Alejandro Ramirez Cajigas]

another simpler graphic example

lat          lon
3.3622      -76.5243
3.3498      -76.5121

z

%points to evaluate
lat	            lon
3.3534876	-76.523249
3.3531877	-76.52344
3.3597285	-76.5240052
3.3591999	-76.5241002

The idea is that I paint and piss within the parameters of the rectangle given by the initial coordinates, currently plotting and measuring within the maximum limit of the area of the outermost circuit, generating the error that it has more area than necessary, because in some cases the boundary of the rectangular wall can cut a circle and for other cases the circles are in the center of the image and external space is lost.

For example now I measure the% are of area adding and subtracting the black and white area, if however it is an error because the black area should be delimited by the coordinates of the walls of the rectangle

            %E = J{:,:}
            lla = [fliplr(table2array(J)) zeros(size(J,1),1)]
            position = lla2flat(lla,[min(lla(:,1)) min(lla(:,2))],90,0)
            % same as before
            %r = ;
            v5= 300
            r = v5 ;
            C = position(:,1:2); % locations (x,y)
            p = 10; % padding
            s = ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
            C = round(C - [min(C(:,1)) min(C(:,2))] + r + p);
            m = false(s);
            m(sub2ind(s,C(:,2),C(:,1))) = 1;
            m22 = bwdist(m)<=r;
            totalarea = sum(m22(:)) % total number of white pixels = area
            m23 = bwdist(m)>=r;
            totalarea2 = sum(m23(:))
            imshow(m22)
            
            
            m24=(totalarea*100)/(totalarea+totalarea2)  %shows area percentage in white
            

you need to put limits with the coordinates of the rectangle so that it would give me the real percentag

                        lat          lon
3.3622      -76.5243
3.3498      -76.5121
%with this he formed the limits of the rectangle
%vertices are (lon1, lat1) (lon2, lat2) (lon2, lat1) (lon1, lat 2)

In example 1 of the stars, there is a reference frame (rectangular), the idea is to measure the area of all the internal stars and the part cut by the reference frame

z

[DGM]

Oof. Sorry about that. Sometimes site notifications get buried, and I almost never check email anymore. Anyway...

% This is only part of the dataset!
lon = [-76.5249196 -76.5281573 -76.5278992 -76.5169571 -76.5322045 -76.5212801 -76.5212291 -76.5090573 -76.5304345 -76.5226856 -76.5180319 -76.5200061 -76.514478 -76.5121552 -76.5076382 -76.5074345 -76.5061631 -76.5059834 -76.5021961 -76.5019253 -76.5001202 -76.4916586 -76.4889863 -76.5249977 -76.5258832 -76.5285806 -76.5167687 -76.5141764 -76.5101736 -76.5053988 -76.4869869 -76.4909459 -76.4945159 -76.5327531 -76.5187478 -76.5101116 -76.5118256 -76.5368692 -76.5358927 -76.5351849 -76.5347021 -76.5364751 -76.5367191 -76.5355765 -76.5315237 -76.4893841 -76.5217632 -76.5215164 -76.5211731 -76.520905 -76.5205857 -76.5203818 -76.5180966 -76.5165489 -76.5165489 -76.5163424 -76.5147331 -76.5147384 -76.5163367 -76.517367 -76.5222486 -76.5227287 -76.523005 -76.5233778 -76.5015155 -76.4964618 -76.4980119 -76.5032478 -76.5232265 -76.5230613 -76.5206634 -76.5264623 -76.525545 -76.5272398 -76.528092 -76.5299112 -76.5320919 -76.5337575 -76.533956 -76.5346507 -76.5356673 -76.5365309 -76.5357638 -76.5348599 -76.5333472 -76.53415 -76.5281195 -76.5268189 -76.5293914 -76.5282836];
lat = [3.4310571 3.4330517 3.4329578 3.4259733 3.4141012 3.4132041 3.4147463 3.4200984 3.4193326 3.4191612 3.4195618 3.4194799 3.4197733 3.4199018 3.4191972 3.4186178 3.4172711 3.4174853 3.4131446 3.413158 3.410971 3.4120227 3.4146176 3.4187462 3.4121893 3.4266255 3.4291025 3.4347968 3.4321028 3.4271533 3.4185984 3.4218167 3.4248689 3.4282086 3.4354939 3.4334187 3.4352097 3.4277076 3.4309797 3.4333837 3.4349393 3.4327974 3.4182723 3.4090406 3.4077261 3.4109503 3.4086604 3.4107276 3.4166607 3.4197919 3.4226876 3.4247333 3.4252902 3.4295177 3.4295177 3.431641 3.434254 3.4347146 3.4324815 3.4213276 3.4175283 3.412915 3.410417 3.4071183 3.4253297 3.4339325 3.4309045 3.4244061 3.4351374 3.4346636 3.4352634 3.4330464 3.4329822 3.432503 3.431776 3.4305833 3.4284118 3.42642 3.4241897 3.4235712 3.4203582 3.4178333 3.4186983 3.4215123 3.4260317 3.431708 3.4312579 3.4321469 3.4104115 3.4129524];
lonlat = table(lon.',lat.');

% boundary
box = [-76.5369 3.4356; -76.4868 3.4071]; % lon/lat

% process table & box; convert to meters
lla = [fliplr(table2array(lonlat)) zeros(size(lonlat,1),1)];
position = lla2flat(lla,[min(lla(:,1)) min(lla(:,2))],90,0);
boxpos = lla2flat([fliplr(box) [0; 0]],[min(lla(:,1)) min(lla(:,2))],90,0);
boxpos = boxpos(:,1:2);

% same as before
r = 400;
C = position(:,1:2); % locations (x,y)

p = 10; % padding
s = ceil([(max(C(:,2))-min(C(:,2))),(max(C(:,1))-min(C(:,1)))] + r*2 + p*2);
refpoint = [min(C(:,1)) min(C(:,2))];
C = round(C - refpoint + r + p);
corners = round(boxpos - refpoint + r + p);

m = false(s);
m(sub2ind(s,C(:,2),C(:,1))) = 1;
m = bwdist(m)<=r;
m(:,[1:corners(1,1) corners(2,1):end]) = 0;
m([1:corners(1,2) corners(2,2):end],:) = 0;

size(C,1)*round(pi*r^2)
totalarea = sum(m(:)) % total number of white pixels = area

imshow(m); hold on
plot(corners([1 2 2 1 1],1),corners([1 1 2 2 1],2),':')

[David Alejandro Ramirez Cajigas]

GRACIAS!!!

Thanks a lot.

I hope one day I can help someone like that. Thanks again