suggestion for geoframe create: read in metadata from .prj file

[rpguiteras]

I have a suggestion for what I think would be a small-but-useful improvement to geoframe create: when loading data from dbf/shp datasets, look for a corresponding .prj file and pull in some metadata.

These could be stored as characteristics in the created .dta file, or in some other way accessible to geoframe and geopoly.

I believe GetSpatialRef() in Python osgeo does something similar:

https://glenbambrick.com/2017/04/24/osgp-data-assessment/

I have added a very simple example below, which could be adapted into an ado-file or a geoframe create option. Unfortunately I have not found any authoritative information on how .prj files are structured, what the possible fields are, etc.

version 18.0
clear all 
cap frame drop nuts03*
cap rm nuts03.dta 
cap rm nuts03_shp.dta 

//# download data
loc GIS_URL "https://github.com/asjadnaqvi/The-Stata-Guide/raw/master/GIS"
loc GIS_FILENAME "NUTS_RG_03M_2021_3035_mainland"

loc FORMATS cpg dbf prj qmd shp shx 
foreach FMT of loc FORMATS {
  copy "`GIS_URL'/`GIS_FILENAME'.`FMT'" "NUTS03.`FMT'" , replace 
}
ls NUTS03*

//# create geoframe
geoframe create nuts03 using NUTS03, current replace 
geoframe save , replace
geoframe describe 
desc, f 
frame dir
frame change nuts03_shp 
desc, f 
save nuts03_shp.dta, replace 
frame change nuts03 

//# read in .prj 
tempname prjfile 
file open `prjfile' using "NUTS03.prj", read text 
file read `prjfile' prjText 
file close `prjfile'
loc prjText = trim(`"`prjText'"')

//# add full text of .prj as a dta char
notes _dta : contents of NUTS03.prj :
notes _dta : `prjText'
char _dta[prj] `prjText'

//# add specific items from prj as dta chars
if strpos(`"`prjText'"',"PROJCS")~=0 {
  char _dta[coord_system_type] "Projected"
  
  // (.*?) for lazy matching -- shortest possible match between (1) and (2)
  // (.*) for greedy matching -- longest possible match
  loc matchProjCode = regexm(`"`prjText'"',`"(PROJCS\[")(.*?)(",)"')
  if `matchProjCode'==1 {
    loc proj_code "`=regexs(2)'" 
    char _dta[projection] `proj_code'
  }
  loc matchProjName = regexm(`"`prjText'"',`"(PROJECTION\[")(.*?)("\],)"')
  if `matchProjName'==1 {
    loc projection_name "`=regexs(2)'" 
    char _dta[projection_name] `projection_name'
  }
}
else if strpos(`"`prjText'"',"GEOGCS")~=0 {
   char _dta[coord_system_type] "Geographic"
}
if strpos(`"`prjText'"',"GEOGCS")~=0 {
  loc matchGCS = regexm(`"`prjText'"',`"(GEOGCS\[")(.*?)(",)"')
  if `matchGCS'==1 {
    loc GCS "`=regexs(2)'" 
    char _dta[GCS] `GCS'
  }
}
loc matchMeter = ustrregexm(`"`prjText'"',`"(UNIT\["Meter",)(.*?)(\])"',0)
if `matchMeter'==1 {
  loc unit_meter "`=ustrregexs(2)'" 
  char _dta[unit_meter] `unit_meter'
}
loc matchSpheroid = regexm(`"`prjText'"',`"(SPHEROID\[")(.*?)(",)"')
if `matchSpheroid'==1 {
  loc GCS_spheroid `"`=regexs(2)'"' 
  char _dta[GCS_spheroid] `GCS_spheroid'
}

//# save and exit 
sleep 1000
save nuts03.dta, replace 
notes _dta 
char li _dta[]
geoframe desc 
frame dir 
exit

[benjann]

Yeah, I know, I've been shying away from this up to now. Obtaining projection info would make it possible to check consistency between data or even make data compatible on the fly.