How to find and apply reprojection parameters?

Question

I have a data set in one projection, and a shapefile in another, and I am trying to combine the two. I am using Python. I apologize for posting a "help me transform coordinates" question on this site - there are so many.

I assume the way to do this conversion is with proj4. I am unable to find/understand the proj4 string needed to do this conversion.

The shapefile.prj contains this:

PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",90.0],UNIT["Meter",1.0]]

Is there a simple way to find the proj4 code for this?

And I have a netCDF file that I can plot correctly with the following Python code:

map = Basemap(llcrnrlon=lon[0,0],
              llcrnrlat=lat[0,0],
              urcrnrlon=lon[560,300],
              urcrnrlat=lat[560,300],
              #rsphere=(6378137.00,6356752.3142),\
              resolution='l',
              area_thresh=1E6,
              projection='lcc',
              #lat_1=90, lon_0=0)
              lat_1=71.0,
              lon_0=-39.0)

I would like to keep the Basemap (python) projection, and convert the shapefile contents to that projection.

Edit

Further research and I found this post: Shapefile PRJ to PostGIS SRID lookup table? with a helpful bit of code. Running that and I find the following:

Proj4 is: +proj=stere +lat_0=90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
EPSG is: None

Edit

I found out I can determine the basemap proj string with:

map.proj4string

I'm now using the advice here http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/ to try to convert from one to the other...

Answer 1

The easiest way to get the PROJ.4 parameters is probably with Fiona.

with fiona.open('shapefile.shp') as source:
    crs = source.crs

The crs will now contain a dictionary mapping of pyproj parameters. To get the PROJ.4 string, you can use fiona.crs.to_string(crs).

You can do this with osr as well, but that's not quite as clean and simple.

Answer 2

I have found a solution.

1) To get proj4 info from a shapefile, I use the code here: https://gis.stackexchange.com/a/7615/609

2) To get the proj4 info for an existing Python BaseMap named 'map': map.proj4string.

3) To set up the transform:

import mpl_toolkits.basemap.pyproj as pyproj
shpProj = pyproj.Proj(<proj string from (1)>)
mapProj = pyproj.Proj(map.proj4string)

4) To transform the shapefile from its projection to lat/lon:

lonlat = array(shpProj(shpX,shpY,inverse=True)).T

5) To convert the lat lon to the Basemap projection:

baseXY = np.array(mapProj(lonlat[:,0],lonlat[:,1])).T

Another implementation of the final transform doesn't involve step (2) or defining the mapProj variable, and just uses the Basemap (here named map) directly:

baseXY = np.array(map.projtran(lonlat[:,0],lonlat[:,1])).T

The final variable, baseXY, is the contents of the shapefile in the coordinate of the Basemap (defined in the original question).

Answer 3

You can feed gdalsrsinfo from GDAL package with your prj file definition, and it returns

PROJ.4 : '+proj=stere +lat_0=90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=
m +no_defs '

as correspondent proj4 string.

Your second file has this definition string:

+proj=lcc +lat_1=71.0 +lon_0=-39.0 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Using ogr2ogr, the command line for projecting the shapefile to the second would be

ogr2ogr -f "ESRI Shapefile" target.shp source.shp -t_srs "+proj=lcc +lat_1=71.0 +lon_0=-39.0 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

You don't have to specify the source srs because it will be read from the prj file.

Answer 4

Besides GDAL/Fiona you can alternatively use PyCRS, a package written in pure Python specifically for defining and converting between and projection formats. So in your case:

import pycrs
wktstring = 'PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",90.0],UNIT["Meter",1.0]]'
projobj = pycrs.parser.from_unknown_wkt(wktstring)
proj4string = projobj.to_proj4()

Or if you want to read it directly from the .prj file:

projobj = pycrs.loader.from_file('path/to/projfile.prj')

It is not nearly as well tested as GDAL though, so it might or might not work on your specific projection string. The documentation has more details on functionality.