Working with Gridcells and Borders

[Prachi Singh]

Hi,

I am trying to create a dataset for a Spatial Regression Discontinuity Design analysis.
For this I use shapefiles from (map-master.zip, subfolder - States): http://projects.datameet.org/maps/states/#download
I split the region of India into gridcells and calculate the distance of these gridcells from the border of treated states (a collection of states).

My code has been pasted below for reference.

I face couple of problems in this exercise (which I list below):

1) I would like to identify gridcells which tend to overlap the border i.e. has part of area lying on either side of the border?

2) I would like to remove certain gridcells from analysis i.e. those which belong to the extreme left corner of the highlighted region (sharing a border with Pakistan and thus that segment of the border only has gridcells from India but not from Pakistan). How can I identify those gridcells?

I have seen @Robert Picard posts on these topics, would be great if you can offer your help regarding this.

Thanks alot for your help,
Prachi

************************************************** ************************************************** *********************

clear

cd "...maps-master/States"

*------------------------------------------------------------------*

*STEP-1: Call the shapefile and save it in stata format
shp2dta using Admin2, data("st_attr.dta") coord("st_coord.dta") ///
genid(stid) gencentroids(cc) replace

*------------------------------------------------------------------*

*STEP-2: Create Gridcells
*Creating ~50*50km gridcells
spgrid using "st_coord.dta", ///
resolution(w0.5) shape(square) ///
cells("India-GridCells.dta") ///
points("India-GridPoints.dta") ///
replace compress dots

*Visualize Gridcells
use "India-GridPoints.dta", clear
spmap using "India-GridCells.dta", id(spgrid_id) ///
polygon(data("st_coord.dta") ///
ocolor(red) osize(thin))

*------------------------------------------------------------------*

*STEP-3: Create a polygon of treated states
use "st_attr.dta", clear
ren stid _ID
ren ST_NM state
mergepoly if inlist(state,"Punjab","Haryana","Delhi") using "st_coord.dta", coor("merged_coord.dta") replace

*------------------------------------------------------------------*

*STEP-4: Create border points
**Source: https://www.statalist.org/forums/for...nt-to-a-border
clear
use "merged_coord.dta", clear

rename _X lon
rename _Y lat

// a new polygon starts with missing values for lat/lon
gen poly = sum(lon == .)
gen id = _n
sort poly id
by poly: assert lon == . if _n == 1
by poly: assert lon[2] == lon[_N]
by poly: assert lat[2] == lat[_N]

sum


// distance between consecutive points
gen double lat2 = lat[_n+1]
gen double lon2 = lon[_n+1]
geodist lat lon lat2 lon2, gen(d) mi

// number of points needed to keep the segments at < 0.1 mi
gen n = 1 + int(d/.1)

expand n if !mi(d)
sort id

// interpolate linearly each intermediary points
by id: gen double lat3 = lat + (_n-1) * (lat2 - lat) / n
by id: gen double lon3 = lon + (_n-1) * (lon2 - lon) / n

// use these new coordinates
drop lat lon
rename lat3 lat
rename lon3 lon
replace id = _n

// double check; ignore missing distances (separate polygons)
keep poly id lat lon
gen double lat2 = lat[_n+1]
gen double lon2 = lon[_n+1]
geodist lat lon lat2 lon2, gen(d) mi
assert round(d,.1) <= .1 if !mi(d)


* Clean-up and save

keep poly id lat lon
save policystates_border_points, replace

*------------------------------------------------------------------*

*STEP-5: Calculate distance of GridCells from the border
use "India-GridPoints.dta", clear
geonear spgrid_id spgrid_ycoord spgrid_xcoord using "policystates_border_points", n(id lat lon) long near(1)
save distance_gridcells, replace

*Visualize gridcells along the border (within 50km)
use "India-GridPoints.dta", clear
merge 1:1 spgrid_id using distance_gridcells
local x 50
gen within`x' = 1 if km_to_id <= `x'
recode within`x'(.=0)
spmap within`x' using "India-GridCells.dta", id(spgrid_id) ///
polygon(data("merged_coord.dta") ///
ocolor(red) osize(thin))

************************************************** ************************************************** *********************

[Scott Merryman]

Prachi Singh I think you can use Robert Picard's -geoinpoly- to find where the borders are within the grid cells.

After your "*STEP-2: Create Gridcells"

clear*
use st_coord
drop _ID
geoinpoly _Y _X using "India-GridCells.dta", noprojection inside
keep if _ID != .
duplicates drop _ID, force
merge 1:m _ID using india-gridcells
keep if _m == 3
drop _m
gen spgrid_id = _ID
duplicates drop
merge 1:m spgrid_id using "India-GridPoints.dta"
keep if _m == 3
spmap using "India-GridCells.dta", id(spgrid_id)  polygon(data("st_coord.dta")  ocolor(red) osize(thin))
save Border_GridCells,replace

[Prachi Singh]

Thanks alot Scott Merryman ! This works perfectly