-1

I am working with historical wind track data, which can be found here.

How can I use pyshp to retrieve all the lat & lon along the windtrack correctly?

I followed the documentation on PyShp and imported 'lin.shp' file type intially but the coordinates returned are not correct when manually checked on Google Earth.
Second, I imported the 'pts.shp' filetype and when I try to run the 'bbox' function it returns an AttributeError: 'Shape' object has no attribute 'bbox'

ah_
  • 1
  • 1

2 Answers2

1

This code should do it

lin_shp = shapefile.Reader("al212021_best_track/AL212021_lin")
coords = [s.points for s in lin_shp.shapes()]
lin_shp.close()

Contents of al212021_best_track dir look like this

 .

 ..

 AL212021_lin.dbf

 AL212021_lin.prj

 AL212021_lin.shp

 AL212021_lin.shp.xml

 AL212021_lin.shx

 AL212021_pts.dbf

 AL212021_pts.prj

 AL212021_pts.shp

 AL212021_pts.shp.xml

 AL212021_pts.shx

 AL212021_radii.dbf

 AL212021_radii.prj

 AL212021_radii.shp

 AL212021_radii.shp.xml

 AL212021_radii.shx

 AL212021_windswath.dbf

 AL212021_windswath.prj

 AL212021_windswath.shp

 AL212021_windswath.shp.xml

 AL212021_windswath.shx

https://pypi.org/project/pyshp/#reading-geometry

AJP
  • 11
  • 3
0

You are retrieving the coordinates correctly. Google Earth requires all data to be EPSG:4326 - WGS84 Geographic. NOAA converts this data for the KML files however the shapefiles are in an esoteric projection which QGIS identifies as "Unknown datum based upon the Authalic Sphere - Projected" for the shapefiles on that page. The map halfway down the data download page you sent is probably in that projection.

Your options are:

  1. Download the KMZ, use python to unzip it, and parse the KML file inside using Python's built-in XML tools to extract the EPSG:4326 points.
  2. Figure out the math to do the transformation of the points yourself to whatever projection you need to work in. Pure Python but this will be difficult.
  3. Switch from using pure Python to using Fiona/Shapely to reproject the points.
  4. Go one level lower from Fiona/Shapely and use the GDAL/OGR Python bindings to reproject the points.
  5. Use GDAL/OGR command line tools called from Python to reproject the points.

Your project and environment restrictions plus your familiarity with these tools will determine which one of those options are best.

Based on the error you are getting on the pts.shp bbox, the problem is you're trying to get a bounding box on a single shape record which contains only one point. A bounding box would only work on shapefile types that have at least two points per record such as the line shapefile in this dataset. You can always get the bounding box at the file level on every shapefile type.

GeospatialPython.com
  • 241
  • 1
  • 6