*! clippolygon v2.0 Asjad Naqvi 08Sep2022

* v1.0: first release. Consolidation of different packages
* Sutherland-Hodgman polygon clipping algorithm
* Asjad Naqvi (asjadnaqvi@gmail.com)



cap program drop clippolygon


********************
// 	clippolygon   //
********************

program define clippolygon, // sortpreserve eclass
	version 15

	syntax namelist(max=1), ///
		Method(string) 	///  				// method is box or circle
		[ Box(numlist min=4 max=4) ] /// 	// box options
		[ XMid(real 0) YMid(real 0) Radius(real 5) Points(real 60) Angle(real 0)  ]  // circle options

	*** checks

	// check gtools
	capture findfile gtools.ado
	if _rc != 0 {
		di as err "gtools package is missing. Click here to install: {stata ssc install gtools, replace}"
		exit
	}

	// carryforward
	capture findfile carryforward.ado
	if _rc != 0 {
		di as err "carryforward package is missing. Click here to install: {stata ssc install carryforward, replace}"
		exit
	}


	capture confirm file "`namelist'.dta"
	if _rc {
	   di as err "{p}File {bf:`namelist'} not found{p_end}"
	   exit 601
	}


	capture confirm file "`namelist'.dta"
	if _rc {
	   di as err "{p}File {bf:`namelist'} not found{p_end}"
	   exit 601
	}

	if ("`method'" != "box" & "`method'" !="circle") {
		di as error "Valid method options are {it:box} or {it:circle}."
	}


	// push the bounds to locals

	if ("`method'" == "box")  {


		if "`box'" == "" {
			di as error "The option {it:box()} must be specified."
			exit 601
		}


	tokenize `box'
		local xmin = `1'
		local xmax = `2'
		local ymin = `3'
		local ymax = `4'

		if `xmax' < `xmin' {
			di as error "xmax < xmin. Please make sure the coordinates are specified in the correct order."
			exit
		}

		if `ymax' < `ymin' {
			di as error "ymax < ymin. Please make sure the coordinates are specified in the correct order."
			exit
		}

		mata: box = `xmin', `ymin' \ `xmin', `ymax' \ `xmax', `ymax' \ `xmax', `ymin' \ `xmin', `ymin'

	}


	if ("`method'" == "circle") {

			mata: box = returnbounds(`xmid', `ymid', `radius', `angle', `points')

			mata: st_local("xmin", strofreal(min(box[.,1])))
			mata: st_local("xmax", strofreal(max(box[.,1])))

			mata: st_local("ymin", strofreal(min(box[.,2])))
			mata: st_local("ymax", strofreal(max(box[.,2])))

	}



// main routine

qui {
	preserve
		use "`namelist'.dta", clear


		// separate out the islands

		cap drop markme
		cap drop group

		sort _ID shape_order
		bysort _ID: gen markme = _n if _X==.

		sort _ID markme
		gen group = _n if markme!=.

		drop markme
		sort _ID shape_order
		carryforward group, replace

		cap drop box
		gen  box = 1 if inrange(_X, `xmin', `xmax') & inrange(_Y, `ymin', `ymax')

		// mark empty _IDs and throw them out

		cap drop temp
		by _ID: egen temp = mean(_X)
		drop if temp==.
		drop temp

		// throw out the shapes that are completely outside the box

		bysort group: egen double boxcheck = mean(box)
		drop if boxcheck==.
		drop boxcheck
		drop box
		sort _ID shape_order

		count

		if `r(N)' == 0 {
			di as error "No shapes are contained within the clipping boundary. Please check the bounds."
			exit
		}


		// mark shapes that are completely inside the box and leave them as they are

		cap drop markme
		gen markme = .


			tempvar minx maxx miny maxy

			bysort group: egen double `minx' = min(_X)
			bysort group: egen double `maxx' = max(_X)

			bysort group: egen double `miny' = min(_Y)
			bysort group: egen double `maxy' = max(_Y)


				if ("`method'" == "box") {
					qui replace markme = 1 if (`minx' >= `xmin') & (`maxx' <= `xmax') &  (`miny' >= `ymin') & (`maxy' <= `ymax')
				}


				if ("`method'" == "circle") {


					local inradius = `radius' * cos(_pi / `points')

					tempvar mark
					gen `mark' = 0

					replace `mark' = `mark' + 1 if (sqrt((`maxx' - `xmid')^2 + (`maxy' - `ymid')^2) < `inradius')
					replace `mark' = `mark' + 1 if (sqrt((`minx' - `xmid')^2 + (`maxy' - `ymid')^2) < `inradius')
					replace `mark' = `mark' + 1 if (sqrt((`maxx' - `xmid')^2 + (`miny' - `ymid')^2) < `inradius')
					replace `mark' = `mark' + 1 if (sqrt((`minx' - `xmid')^2 + (`miny' - `ymid')^2) < `inradius')

					qui replace markme = 1 if `mark'==4

				}



		// save here as a tempfile since we cannot preserve twice

		drop `minx' `maxx' `miny' `maxy'
		cap drop `mark'

		sort _ID shape_order group
		tempfile _polysample
		save "`_polysample'", replace


		// store the shapes that are fully inside in a separate file. This speeds up teh calculations
		keep if markme==1
		drop markme
		tempfile _polyinside
		save "`_polyinside'", replace


		// continue with the shapes that intersect with the box
		use "`_polysample'", clear


		di "Values to be dropped"
		count if markme==1


		drop if markme==1
		drop markme


		di "Values left"
		count


		// start the procedure here. Each island needs to be processed separately.

		cap gen double clip_x = .
		cap gen double clip_y = .
		cap gen double id = .
		cap gen double group = .


		mata: points   = st_data(., ("_X", "_Y", "_ID", "group"))
		mata: myindex  = uniqrows(points[.,4])
		mata: myindex  = select(myindex, (myindex[.,1] :< .))
		mata: clipbox  = returnclip(points, myindex, box)
		mata: st_local("newobs", strofreal(rows(clipbox)))


		// expand observations (automate this for the custom data range)
			if `newobs' > _N {
				set obs `newobs'
			}

			if `newobs' > 1 {
				getmata (clip_x clip_y id group shape_order) = clipbox, force double replace
			}


		keep clip_x clip_y id group shape_order

		drop if id==.
		ren id _ID
		ren clip_x _X
		ren clip_y _Y
		gen clipped = 1

		append using "`_polyinside'" // add back the inside polygons



		order _ID _X _Y group shape_order
		sort _ID group shape_order

		compress
		save "`namelist'_clipped.dta", replace
	restore
 }

	di in green "Done! File exported as `namelist'_clipped.dta"


end


**************************
//   Mata subroutines   //
**************************


******************
//  returnclip  //
*****************

cap mata: mata drop returnclip()

mata:
function returnclip(data, index, box)
{
	mycoords = J(1,5,.)

	for (i=1; i <= rows(index); i++) {


		mysubset = select(data, data[.,4] :== index[i] )
		mysubset = select(mysubset, (mysubset[.,1] :< .))

		myclips = clipme(mysubset, box)

		if (rows(myclips) > 1) {

			myclips[1,3] = myclips[2,3]
			myclips[1,4] = myclips[2,4]
			myclips =  myclips   \ ., .,  myclips[2,3], myclips[2,4] //   ., .,  myclips[2,3], myclips[2,4]

			myclips = myclips, (1::rows(myclips))
			mycoords = mycoords \ myclips

		}

	}

	return (mycoords)
}
end


**********************
//  returnbounds	//
**********************

cap mata: mata drop returnbounds()
mata:  // returnbounds
	function returnbounds(x, y, r, a, o)  // xmid, ymin, radius, angle, obs
	{
		theta  = J(o,1,.)

		for (i=1; i <= o; i++) {
			theta[i] = i * -2 * pi() / o
		}

		coords = cos(theta) :* r, sin(theta) :* r
		ro = -1 * a * pi() / 180
		rotation = (cos(ro), -sin(ro) \ sin(ro), cos(ro))
		coords = (coords * rotation') :+ (x, y)
		coords = coords \ coords[1,.]

		return (coords)

	}
end

********************
//  is_inside	  //
********************

cap mata: mata drop is_inside()
mata:  // is_inside
	function is_inside(p1,p2,q)  // clip point start, clip point end, point to be evaluated
	{
		scalar R
		R = (p2[.,1] - p1[.,1]) * (q[.,2] - p1[.,2]) - (p2[.,2] - p1[.,2]) * (q[.,1] - p1[.,1])
		return (R <= 0)
	}

end


****************************
//  compute_intersection  //
****************************

cap mata: mata drop compute_intersection()
mata:  // compute_intersection

	function compute_intersection(p1,p2,p3,p4)  //
	{

		scalar x, m2, b2, y, m1, b1
		vector intersection

        // if first line is vertical
		if (p2[.,1] - p1[.,1] == 0) {

			x = p1[.,1]
			m2 = (p4[.,2] - p3[.,2]) / (p4[.,1] - p3[.,1])
			b2 = p3[.,2] - m2 * p3[.,1]
			y = m2 * x + b2
        }

        // if second line is vertical
		else if (p4[.,1] - p3[.,1] == 0) {

			x = p3[.,1]
			m1 = (p2[.,2] - p1[.,2]) / (p2[.,1] - p1[.,1])
			b1 = p1[.,2] - m1 * p1[.,1]
			y = m1 * x + b1
        }

        // if neither line is vertical
        else {

			m1 = (p2[.,2] - p1[.,2]) / (p2[.,1] - p1[.,1])
			b1 = p1[.,2] - m1 * p1[.,1]

			m2 = (p4[.,2] - p3[.,2]) / (p4[.,1] - p3[.,1])
			b2 = p3[.,2] - m2 * p3[.,1]

			x = (b2 - b1) / (m1 - m2)
			y = m1 * x + b1
		}

		intersection = (x, y, p1[.,3], p1[.,4])
        return (intersection)
	}
end



**************
//  clipme  //
**************

cap mata: mata drop clipme()
mata:  // clipme

	function clipme(subpoly,clippoly)
	{

	scalar i, j, cedge1, cedge2, sedge1, sedge2
	matrix finalpoly, nextpoly, intersection

	finalpoly = subpoly


		for (i=1; i <= rows(clippoly) - 1; i++) {

			if (finalpoly[.,1]!=.) {
				finalpoly   = select(finalpoly, (finalpoly[.,1] :< .))
			}

			finalpoly   = finalpoly \ finalpoly[1.,]
			nextpoly = finalpoly
			finalpoly = J(1, 4, .)

			cedge1 = clippoly[i    , .]
			cedge2 = clippoly[i + 1, .]

			for (j=1; j <= rows(nextpoly) - 1; j++) {

				sedge1 = nextpoly[j    , .]
				sedge2 = nextpoly[j + 1, .]

				if (is_inside(cedge1,cedge2,sedge2) == 1) {
					if (is_inside(cedge1,cedge2,sedge1) == 0) {
						intersection = compute_intersection(sedge1, sedge2, cedge1, cedge2)
						finalpoly = finalpoly \ intersection
					}
					finalpoly = finalpoly \ sedge2
				}
				else if (is_inside(cedge1,cedge2,sedge1) == 1) {
					intersection = compute_intersection(sedge1, sedge2, cedge1, cedge2)
					finalpoly = finalpoly \ intersection
				}
			}
		}

		// if have a line, clean it up
		if (rows(finalpoly) == 4) {
			if (finalpoly[2,1]!=finalpoly[3,1] & finalpoly[2,2]!=finalpoly[3,2]) {
				return (finalpoly[1::rows(finalpoly) - 1, .])
			}
		}
		else {
			return (finalpoly)
		}
	}
end

****************************
***     END OF FILE      ***
****************************