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...

score 3 · Answer 1 · answered Nov 14 '13 at 18:09

3

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.

answered Nov 14 '13 at 18:09

Kelsey

331
1
3

So the crs string from this code is '+k=1 +lat_0=90 +lon_0 +no_defs +x_0 +y_0' which is quite different than the proj4 string I found using the osr and that function provided here http://gis.stackexchange.com/questions/7608/shapefile-prj-to-postgis-srid-lookup-table – mankoff Nov 14 '13 at 18:28
That definition looks a bit broken, with no ellipsoid and no projection method declared. – AndreJ Nov 14 '13 at 19:10

score 3 · Accepted Answer · edited Apr 13 '17 at 12:33

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).

AndreJ · Answer 3 · 2013-11-14T20:24:34.840

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.

Thanks. That appears to be the projection I found (see Edit). I am now working on converting that to the basemap projection (see 2nd edit)... — mankoff, Nov 14 '13 at 19:18

score 0 · Answer 4 · answered Aug 29 '15 at 13:37

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.

How to find and apply reprojection parameters?

Edit

Edit

4 Answers4