3

I am trying to work out how to apply WGS84 to a POLYGON I am creating from a list of latitude and longitude coordinates. Without the correct coordinate system, I cannot calculate the area, my intersection is also failing.

# Create Polygon from formatted string of coordintes
wkt1 = "POLYGON (({0}))".format(coords_fin)
poly1 = ogr.CreateGeometryFromWkt(wkt1)

I know I can define the WGS84 geometry and have successfully done this for other tasks:

# Set up the coordinate reference system, WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

What I can't find is how to say my polygon is in this CRS. All the examples I have managed to find are based on exporting results as shapefiles, etc, so are applying the CRS to a layer.

I simply want to define two polygons in WGS84 (sanity check the area), then return a False if they don't overlap using:

intersection = poly1.Intersection(poly2)

I have tried SetGeometry but it doesn't appear to be an attribute of the polygon.

1
  • But does the intersection return what you expect? If they are both in the same SRS then the intersection should work. You only need the SRS specification if you need to transform either layer. Dec 20 '16 at 10:44
2

The wonderful Python GDAL/OGR Cookbook has a whole section on "Projection": https://pcjericks.github.io/py-gdalogr-cookbook/projection.html

The main points are these:

You need have both a source and target system:

source = osr.SpatialReference()
source.ImportFromEPSG(2927)
target = osr.SpatialReference()
target.ImportFromEPSG(4326)

Then you can create a Transformation between those systems:

transform = osr.CoordinateTransformation(source, target)

And finally apply the transformation to geometries:

point.Transform(transform)
5
  • But the list of coordinates I'm using are in WGS84 already. I am not trying to change the projection, just tell gdal to interpret my polygon correctly.
    – emmalg
    Dec 16 '16 at 14:22
  • In other words both the source and the destination would have the same projection.
    – emmalg
    Dec 16 '16 at 17:45
  • Oops- pressed enter... The list of source coordinates were recorded in WGS84 and that's what I need to apply. I tried applying a WGS84-WGS84 transformation and although it didn't complain it didn't seem to make any difference..
    – emmalg
    Dec 16 '16 at 17:47
  • How does GDAL not interpret your polygon correctly? Dec 18 '16 at 21:47
  • Both polygons are in Greenland and one is entirely within the other. I was not getting the expected results using GDAL when checking the intersection (or I was trying to do something wrong). However, it is a moot point now, I urgently had to clear TBs of space on a disk before I went on holiday on Friday. I decided to try Basemap which I'm more familiar with. I set up WGS84 before loading the data. I thought it would be too slow. Lucky it was really quick to do the intersection and I was able to get rid of the unwanted data easily.
    – emmalg
    Dec 21 '16 at 22:03

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.