Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • New on SSC: Create circles defined by geographic coordinates with -geocircles-

    Thanks to Kit Baum, a new command called geocircles is now available from SSC.

    geocircles creates polygons that approximate the shape of circles drawn on a sphere. Each polygon point is at a great-circle distance from the center point, defined by its latitude and longitude. A great-circle distance is the shortest distance between two points measured along the surface of a sphere.

    geocircles generates two datasets that follow the model used by shp2dta (from SSC) when converting shapefiles to Stata datasets.

    Circles generated by geocircles will appear deformed if plotted using unprojected coordinates. Use geo2xy (from SSC, see announcement here) to apply a map projection to geographic coordinates before creating a map in Stata.

    Thanks to Luca Aguzzoni who asked a question a few weeks ago about creating a map that would show a point and a circle around it. My proposed solution evolved into geocircles.

    To install geocircles, type in Stata's command window

    Code:
    ssc install geocircles

  • #2
    Hello Mr. Picard,
    I am trying to use the twoway graph commands to produce a map on which I can draw a series of circles of 10 miles from a series of points (geographic cordinates for a railroad). geocircles seems like it would be the perfect program to complete this task. However, my data is in USA Contiguous Albers Equal Area Conic projection and has the following projection:

    PROJCS["USA_Contiguous_Albers_Equal_Area_Conic",GEOGC S["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-96.0],PARAMETER["Standard_Parallel_1",29.5],PARAMETER["Standard_Parallel_2",45.5],PARAMETER["Latitude_Of_Origin",37.5],UNIT["Meter",1.0]]

    Here is an example of what the data looks like:
    _ID xcord ycord
    1
    1 1703095.1 643192.19
    1 1703165.3 643179.33
    1 1703223.2 643139.69
    1 1703227.6 643134.58
    2
    2 1703227.6 643134.58
    2 1703271.1 643086.92
    2 1703596.4 642574.43
    2 1703677.6 642464.74
    2 1703785.9 642291.48
    3
    3 1644111.8 371236.96
    3 1644166.1 371260.89
    4
    4 1644111.8 371236.96
    4 1643982.1 371365.1

    Geocircles appears to require the data to be latitude and longitude. Are either of the three following solutions possible:
    1. (Preferred method) Convert from the USA Contiguous Albers Equal Area Conic projection to latitude and longitude. I know it is possible to go from lat/lon to my data projection using the following command: geo2xy ycord xcord , replace proj(albers, 6378137 298.257222101 29.5 45.5 37.5 -96)
    2. Run geocircles with the USA Contiguous Albers Equal Area Conic projection instead of latitude and longitude.
    3. Modify the program you posted here: https://www.statalist.org/forums/for...les-around-dot to use the equal albers equal area conic projection.

    Comment


    • #3
      Your best bet is to convert your data to lat/lon and then use geocircles (from SSC) to create circle polygons. Unfortunately, geo2xy (from SSC) does not convert projected coordinates to lat/lon. You would not ask this question if you had access and experience with ArcGIS so I assume that's not an easy option for you. I do not have any experience with other GIS programs but perhaps others can chime in.

      In the mean time, you may want to take a look at a batch converter from mygeodata. Click on the "Choose input coordinate systems" button and type "Albers". Select the "5070 NAD83/ Conus Albers" option (should be the third in the list). You need to change the lat_0 parameter to match your projection. The field should read:

      HTML Code:
      +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
      I then pasted the following input coordinates (the first 4 points from your data example):
      Code:
      1703095.1 643192.19
      1703165.3 643179.33
      1703223.2 643139.69
      1703227.6 643134.58
      For the Output coordinate system I picked the first one on the list ("4326 WGS 84").The field shows:
      Code:
      +proj=longlat +datum=WGS84 +no_defs
      After I clicked on the "Transform" button, I got the following results:
      Code:
      -75.2585854772;41.5752114386
      -75.2577910738;41.5749634344
      -75.2572118236;41.5745056123
      -75.2571732929;41.5744525225
      I copied these lat/lon coordinates to Stata and used geo2xy to convert them using the same projection. The results match what I started from exactly:

      Code:
      . geo2xy lat lon , gen(y x) proj(albers, 6378137 298.257222101 29.5 45.5 37.5 -96)
      
      . list lat lon x y
      
           +------------------------------------------------+
           |       lat          lon           x           y |
           |------------------------------------------------|
        1. | 41.575211   -75.258585   1703095.1   643192.19 |
        2. | 41.574963   -75.257791   1703165.3   643179.33 |
        3. | 41.574506   -75.257212   1703223.2   643139.69 |
        4. | 41.574453   -75.257173   1703227.6   643134.58 |
           +------------------------------------------------+

      Comment


      • #4
        I am not a expert user in stata but working on a coordinates dataset for district boundaries in Pakistan and trying to create x and y cartesian coordinates for point/labels in those districts. Purpose of this exercise to develop choropleth maps of districts and labels on districts. I can create choropleth maps but unable to label districts because my dataset dont have x and y coordinates for district centers.



        use "uscoord.dta", clear
        geo2xy _Y _X, projection( mercator) replace
        local lon0 `r(lon0)'
        local f `r(f)'
        local a `r(a)'
        save "uscoord_mercator.dta", replace

        use usdb,clear
        gen x = runiform()
        save "usdb2",replace

        gen labtype = 1
        append using usdb2
        replace labtype = 2 if labtype==.
        replace ADM2_EN = string(_ID, "%3.0f") if labtype ==2
        gen x_c=runiform()
        gen y_c=runiform()

        keep x_c y_c ADM2_EN labtype
        geo2xy y_c x_c, projection( mercator, `a' `f' `lon0' ) replace
        save maplabels, replace

        clear
        import excel "corana_cases.xlsx", sheet("District") firstrow
        merge 1:1 _ID using usdb2
        drop if _merge!=3

        spmap cases using uscoord, id(_ID) fcolor(BuRd) clmethod(custom) clbreaks(1 10 20 30 40 500) ndlabel("No Cases") legend(position(5)) title("District Wise Confirmed COVID-19 Cases") subtitle("26 March 2020") ///
        label(data(maplabels) xcoord(x_c) ycoord(y_c) ///
        label(ADM2_EN) by(labtype) size(*0.85 ..) pos(12 0) )

        But when i run this code, x and y coordinates are generated in negative values and map is not created.

        I will be thankful for support.

        Comment

        Working...
        X