0

I am a newbie working on GIS data.

I am working on the problem of finding overlapping regions using polygon shapefiles. Most of my shapefiles contain information with latlongs. For example the following is an example of correct polygon shapefile.

'POLYGON ((-122.2787982036204 37.82711080504187, -122.2785692208205 37.82702422544803, -122.2784264333931 37.82771308306625, -122.2786006249591 37.82834442187075, -122.278618196814 37.82840715086107, -122.2786275307938 37.82844027911263, -122.2786494440236 37.8285180610707, -122.2786803950343 37.82862874157701, -122.2787505712275 37.8288838276107, -122.2785544240381 37.82944163165736, -122.2784231200261 37.82977186363865, -122.2782261706852 37.83026719288048, -122.2779896930518 37.83092718669832, -122.2778674662716 37.83126830979819, -122.2777110030613 37.83124252122994, ... )

However, one of the shap file do not contain latlong but instead looks like this:

POLYGON (((-13623529.36072154 6065172.989204641, -13622963.59978288 6065082.239418259, -13622939.17208787 6065083.360996598, -13622930.12814104 6065083.681784724, -13622921.08444332 6065083.807151197, -13622909.33942664 6065083.65457371, -13622902.09806475 6065083.431633831, -13622894.40669532 6065083.099930355, -13622877.68339248 6065082.100505441, -13622874.09881658 6065081.641047121, -13622867.41827274 6065080.370038355, -13622864.73663 6065079.939777289, -13622858.42273742 6065079.249990978, -13622852.1025931 6065078.774796806, -13622846.22608016 6065078.527636741, -13622840.34040108, ...)))

I searched a lot but I could not find what these numbers represent. Is there a way to convert this shapefile into the other format with lat-longs. If not, then how to find the overlapping area of this shapefile with other shapefile with latlongs.

I am currently using the following code in python using shaply:

layer1 = row1["geometry"] # shapefile1
layer2 = row2["geometry"] # shapefile2
layer3 = layer1.intersection(layer2) # overlapping layer

This code is working well for all shapefiles in the latlong format but not with the one have these long numbers.

  • 1
    Welcome to GIS StackExchange! It looks like your second set of data is in a different coordinate system/projection. Where did you obtain the data? They should have some metadata stating what the projection is, otherwise you may need to experiment to identify the current projection, and then reproject into the coordinate system that matches the rest of your data. – smiller Nov 14 '18 at 16:50
  • Thanks smiller, I followed from Ian's answer and identified the CRS of the data. Now its working fine. – userxxx Nov 14 '18 at 17:03
4

The shapefile with the long numbers in is in a different projection (CRS, SRS, Coordinate System) than the others so it's units are probably metres (or feet or US Survey feet or chains, or some such). Until you find out what it's projection is there is nothing you can do with the data as it is just a bunch of (somewhat correlated) random numbers.

So, look in your data directory and see if there is a .prj file with the same basename as your shapefile. This will tell you (and more importantly shapely) what the numbers mean.

Then you can do something like (from the docs):

from functools import partial
import pyproj
from shapely.ops import transform

project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:4326'), # source coordinate system
    pyproj.Proj(init='epsg:26913')) # destination coordinate system

g2 = transform(project, g1)

And if you don't know the EPSG code then you can go directly from the .prj file according to this question using:

prjText = open( prjPath, 'r').read()
srs = osr.SpatialReference()
if ( srs.ImportFromWkt( prjText ) ):
    return ( False, "error importing .prj information from " + prjPath )
inProjection = pyproj.Proj( srs.ExportToProj4() )
  • Thank you so much @Ian Turton, I checked the meta data and identified the CRS. and did this in my code : state_cts = state_cts.to_crs(epsg='4326'). It is working perfectly fine now. Thanks :) – userxxx Nov 14 '18 at 17:02

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.