National Oceanic and
Atmospheric Administration
United States Department of Commerce

Using Shapefiles in PyFerret

The U.S. government is closed. This site will not be updated; however, NOAA websites and social media channels necessary to protect lives and property will be maintained. To learn more, visit commerce.gov.

For the latest forecasts and critical weather information, visit weather.gov.

Using Shapefiles in PyFerret

Question:

Can I use Shapefiles in Ferret or PyFerret?

 

Solution:

PyFerret has a set of SHAPEFILE* functions to read and write shapefiles.  (Try SHOW FUNCTION or SHOW FUNCTION/DETAILS for more information).

yes? show function/brief shape*

SHAPEFILE_READXY(SHAPEFILE,MAXPTS)
SHAPEFILE_READXYVAL(SHAPEFILE,VALNAME,MAXPTS)
SHAPEFILE_READXYZ(SHAPEFILE,MAXPTS)
SHAPEFILE_READXYZVAL(SHAPEFILE,VALNAME,MAXPTS)
SHAPEFILE_WRITEVAL(SHAPEFILE,VALUE,VALNAME,MAPPRJ)
SHAPEFILE_WRITEXYVAL(SHAPEFILE,GRIDX,GRIDY,VALUE,VALNAME,MAPPRJ)
SHAPEFILE_WRITEXYZVAL(SHAPEFILE,GRIDX,GRIDY,GRIDZ,VALUE,VALNAME,MAPPRJ)

These use the PyShp Python package, so you will need to make sure that PyShp is installed.  The PyShp package (see: https://pypi.org/project/pyshp/ ) is pure Python code and can be installed using "pip" or possibly through your system's package manager (We have seen it under the package names pyshp and python-pyshp).

 

Example:

A basic example to get us started, showing use of the SHAPEFILE_READXY function.  Download the 2010 US Counties shapefile information, which consists of the following files:
 
    tl_2010_us_county10.dbf
    tl_2010_us_county10.prj
    tl_2010_us_county10.shp
    tl_2010_us_county10.shp.xml
    tl_2010_us_county10.shx
 

Read and work with this shapefile using these commands from PyFerret's "ferret" command line:

yes? set mem /size=500
yes? let uscounties = shapefile_readxy("tl_2010_us_county10", 7560000)

yes? load /perm uscounties
yes? let uscounties_poslon = if uscounties ge 0.0 then uscounties else 360.0 + uscounties
yes? plot /line=7 /vs /hlimits=235:250 /vlimits=42:50 uscounties_poslon[j=1],uscounties[j=2]
Here we happened to know there are no more that 7560000 coordinates, so using that for argument 2 of SHAPEFILE_READXY lets the function call return immediately, without reading anything until the LOAD is done.  If you did not know the count, you can use -1 and it will do an initial read of the shapefile to get the number of coordinates (so it will take a while to do the initial read). 
 
You may need to increase the upper bound of allowed memory, as was done here, if you are reading a large shapefile.  Note that shapefile coordinates are just numbers, not necessarily longitude and latitude.  The .prj file is there to explain these numbers, but we have not delved into know how much the PyShp package makes use of this file.  The returned value is an array of X,Y coordinates, where X values are at [J=1] and Y values are at [J=2].  The plot in this example is just an ordinary X-Y plot (after adjusting the X values - which happen to be longitudes in this case - to positive values).