Cylindrical Equidistant WRF projection not supported

[SHEN-Cheng]

Hi,
When i use shapefile to subset a wrfout data, and i got an error

   177         raise NotImplementedError('WRF proj not implemented yet: 
--> 178                                   '{}'.format(proj_id))
    179 
    180     proj = gis.check_crs(p4)

NotImplementedError: WRF proj not implemented yet: 6

I checked my data , the projection is:

MAP_PROJ :6     
MAP_PROJ_CHAR :    Cylindrical Equidistant

Is there any way to solve this problem, this projection is the most common projection for NetCDF data

[fmaussion]

Thanks for the report.

this projection is the most common projection for NetCDF data

Yes, but definitely not for WRF data. I actually hear of someone using this projection with WRF for the first time.

Is there any way to solve this problem

Yes. We accept pull-requests. The change probably needs to be implemented here:

salem/salem/sio.py

Lines 146 to 177 in 0a0a930

	# Normal WRF file
	cen_lon = ds.CEN_LON
	cen_lat = ds.CEN_LAT
	dx = ds.DX
	dy = ds.DY
	pargs['lat_1'] = ds.TRUELAT1
	pargs['lat_2'] = ds.TRUELAT2
	pargs['lat_0'] = ds.MOAD_CEN_LAT
	pargs['lon_0'] = ds.STAND_LON
	pargs['center_lon'] = ds.CEN_LON
	proj_id = ds.MAP_PROJ
	if proj_id == 1:
	# Lambert
	p4 = '+proj=lcc +lat_1={lat_1} +lat_2={lat_2} ' \
	'+lat_0={lat_0} +lon_0={lon_0} ' \
	'+x_0=0 +y_0=0 +a=6370000 +b=6370000'
	p4 = p4.format(**pargs)
	elif proj_id == 2:
	# Polar stereo
	p4 = '+proj=stere +lat_ts={lat_1} +lon_0={lon_0} +lat_0=90.0' \
	'+x_0=0 +y_0=0 +a=6370000 +b=6370000'
	p4 = p4.format(**pargs)
	elif proj_id == 3:
	# Mercator
	p4 = '+proj=merc +lat_ts={lat_1} ' \
	'+lon_0={center_lon} ' \
	'+x_0=0 +y_0=0 +a=6370000 +b=6370000'
	p4 = p4.format(**pargs)
	else:
	raise NotImplementedError('WRF proj not implemented yet: '
	'{}'.format(proj_id))

It's probably relatively simple to create the PROJ params for this projection out of the file. If you open a PR, please also create a small size version of the file to parse so that we can add it to the salem test files for the future. Thanks!

[SHEN-Cheng]

Thanks for the report.

this projection is the most common projection for NetCDF data

Yes, but definitely not for WRF data. I actually hear of someone using this projection with WRF for the first time.

Is there any way to solve this problem

Yes. We accept pull-requests. The change probably needs to be implemented here:

salem/salem/sio.py

Lines 146 to 177 in 0a0a930

	# Normal WRF file
	cen_lon = ds.CEN_LON
	cen_lat = ds.CEN_LAT
	dx = ds.DX
	dy = ds.DY
	pargs['lat_1'] = ds.TRUELAT1
	pargs['lat_2'] = ds.TRUELAT2
	pargs['lat_0'] = ds.MOAD_CEN_LAT
	pargs['lon_0'] = ds.STAND_LON
	pargs['center_lon'] = ds.CEN_LON
	proj_id = ds.MAP_PROJ
	if proj_id == 1:
	# Lambert
	p4 = '+proj=lcc +lat_1={lat_1} +lat_2={lat_2} ' \
	'+lat_0={lat_0} +lon_0={lon_0} ' \
	'+x_0=0 +y_0=0 +a=6370000 +b=6370000'
	p4 = p4.format(**pargs)
	elif proj_id == 2:
	# Polar stereo
	p4 = '+proj=stere +lat_ts={lat_1} +lon_0={lon_0} +lat_0=90.0' \
	'+x_0=0 +y_0=0 +a=6370000 +b=6370000'
	p4 = p4.format(**pargs)
	elif proj_id == 3:
	# Mercator
	p4 = '+proj=merc +lat_ts={lat_1} ' \
	'+lon_0={center_lon} ' \
	'+x_0=0 +y_0=0 +a=6370000 +b=6370000'
	p4 = p4.format(**pargs)
	else:
	raise NotImplementedError('WRF proj not implemented yet: '
	'{}'.format(proj_id))

It's probably relatively simple to create the PROJ params for this projection out of the file. If you open a PR, please also create a small size version of the file to parse so that we can add it to the salem test files for the future. Thanks!

I referred to wrf-python and pyproj source code and got the result that this projection standard name is"Equidistant Cylindrical (Plate Caree)"("eqc" for short)

    def _proj4(self):
        _proj4 = ("+proj=eqc +units=m +a={} +b={} "
                  "+lon_0={} +nadgrids=@null".format(
                      Constants.WRF_EARTH_RADIUS,
                      Constants.WRF_EARTH_RADIUS,
                      self.stand_lon))
        return _proj4

The above code is the about that projection in WRF-Python.
You can try the code my type, but i do not know the result is true or not.

    elif proj_id == 6:
        # Equidistant Cylindrical (Plate Caree)
        p4 = '+proj=eqc +lat_ts={lat_1} ' \
             '+lon_0={center_lon} ' \
             '+x_0=0 +y_0=0 +a=6370000 +b=6370000'
        p4 = p4.format(**pargs)

I cannot convert it to usable salem code, sorry

[fmaussion]

Thanks for the details! Looks quite right to me yes.

I cannot convert it to usable salem code, sorry

OK, then we'll see if someone else wants to pick that up and solve this issue!

I think it might help future potential contributors if you could provide a minimal netcdf file here, with only the lons lats, eastings northings and netcdf file attributes for testing (example: https://github.com/fmaussion/salem-sample-data/blob/master/salem-test/grid/geo_em_d01_lambert.nc )

[SHEN-Cheng]

Thanks for the details! Looks quite right to me yes.

I cannot convert it to usable salem code, sorry

OK, then we'll see if someone else wants to pick that up and solve this issue!

I think it might help future potential contributors if you could provide a minimal netcdf file here, with only the lons lats, eastings northings and netcdf file attributes for testing (example: https://github.com/fmaussion/salem-sample-data/blob/master/salem-test/grid/geo_em_d01_lambert.nc )

[SHEN-Cheng]

Thanks for the details! Looks quite right to me yes.

I cannot convert it to usable salem code, sorry

OK, then we'll see if someone else wants to pick that up and solve this issue!
I think it might help future potential contributors if you could provide a minimal netcdf file here, with only the lons lats, eastings northings and netcdf file attributes for testing (example: https://github.com/fmaussion/salem-sample-data/blob/master/salem-test/grid/geo_em_d01_lambert.nc )
Sure！
test.zip

[matthiasdemuzere]

I am looking into this for the integration of LCZ classes in WRF's geo_em files. So far I get a mismatch between the domains for the Cylindrical Equidistant projection. If I can get it right, I'll see if I can push the solution to salem as well ...

[matthiasdemuzere]

@fmaussion I wanted to push a potential solution for the MAP_PROJ = 6 case to a new branch, but I am not allowed?
What did I miss?

[fmaussion]

Thanks for the contribution! Have you made a pull-request before? The idea would be to fork salem and then open a PR to the main (master) branch: https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/about-pull-requests