0

I've never worked with shapefiles before and I'm having a lot of trouble with something that I think should be pretty simple.

So I have the census tract files from NHGIS for 1980 and the structure looks like:

NHGISST     NHGISCTY    GISJOIN     GISJOIN2    SHAPE_AREA  SHAPE_LEN   geometry
0   560     0210    G56002100016    56002100016     2.569850e+09    241261.569958   POLYGON ((-732462.5590209048 499607.0109478103...
1   560     0210    G56002100013    56002100013     9.690180e+06    14271.148305    POLYGON ((-730612.5504814894 448835.8542251377...
2   560     0210    G56002100014    56002100014     1.686311e+07    19833.251952    POLYGON ((-728941.3490845333 443959.3627771207...
3   560     0210    G56002100008    56002100008     4.748756e+06    10873.293718    POLYGON ((-733509.072211165 444652.6344217169,...
4   560     0210    G56002100009    56002100009     3.887032e+06    10440.687964    POLYGON ((-734860.1103011537 445048.4184498419...

I am trying to use the geometry polygon structure to extract a lat/lon centroid but nothing seems to work. Whenever I use something like MultiPoint or Proj (I'm using python by the way) I get coordinates that are completely off.

I'm completely stumped.

def to_lat_long(t):
    inProj = Proj(init='epsg:3857')
    outProj = Proj(init='epsg:4326')
    return transform(inProj,outProj,t[0],t[1])
break_count, break_count_max = 0, 1
for k,v in shape_file_store_1.items():
    if type(v[1990]) == str:
        x = v[1990][1:]
        for s in shp_1990:
            if s == 0:                
                first = next(iter(shp_1990))
            else:
                first = s

            shp_file = first['properties']
            items = list(shp_file.items())
            s_file = items[3][1]

            coords = first['geometry']['coordinates'][0]

            points = MultiPoint(coords)
            points_c = points.centroid
            points_r = points.representative_point()

            print(f'Centroid: {points_c}')

            centroid_latlong = to_lat_long(points_c.coords[0])
            df_townhalls.at[df_townhalls['geo'] == v['geo'], '1990_lat'] = centroid_latlong[0]
            df_townhalls.at[df_townhalls['geo'] == v['geo'], '1990_long'] = centroid_latlong[1]

            print(v['geo'])
            #print(f'Representative Point: {points_r}')
            #print(df_townhalls.loc[df_townhalls['geo'] == v['geo']])

            break
            if s_file == x:
                v[1990] = first
        #break
        break_count += 1
        if break_count > break_count_max: break
    #break
  • 1
    code added @BERA – jdv12 Feb 19 '19 at 17:27
  • What output do you get and what do you expect? – BERA Feb 19 '19 at 17:29
  • I am getting lat/lon coordinates that are nowhere near the census tract/not even in the use. I'm looking to lat/lon centroid of each of these tracks – jdv12 Feb 19 '19 at 17:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.