Missing geographic crs from projcrs for MODIS sinusoidal projection

I need to write raster data projected to the MODIS sinusoidal grid to a netCDF CF file. I am able to get a custom projcrs object from a colleague, who used a custom wkt string that describes the sinusoidal grid using projcrs. However, the projcrs.GeographicCRS property is empty:

wkt = "PROJCS[""MODIS Sinusoidal"",BASEGEOGCRS[""User"",DATUM[""World Geodetic Survey 1984"",SPHEROID[""Authalic_Spheroid"",6371007.181,0.0]],PRIMEM[""Greenwich"",0.0],UNIT[""Degree"",0.0174532925199433]],PROJECTION[""Sinusoidal""],PARAMETER[""False_Easting"",0.0],PARAMETER[""False_Northing"",0.0],PARAMETER[""Central_Meridian"",0.0],UNIT[""Meter"",1.0]]";
proj = projcrs(wkt);
proj = 

projcrs with properties:

Name: "MODIS Sinusoidal"
GeographicCRS: [0×0 geocrs]
ProjectionMethod: "Sinusoidal"
LengthUnit: "meter"
ProjectionParameters: [1×1 map.crs.ProjectionParameters]

I need to get the Earth radius of the Authalic sphere, and I'd expect to get that from the Geographic CRS.

I can use the projcrs.wktstring to get the wkt back and then parse that, but it's a pain. Shouldn't the projcrs object contain the GEOGCRS information from the original wkt?

Thanks for any advice.

MJ

Answer 1

[Karan Singh]

Hi Mary, in the case you described, it seems that the “projcrs” object you obtained from your colleague does not populate the “GeographicCRS” property. This property should contain the Geographic CRS information from the original “WKT” (Well-Known Text) string.

The absence of the “GeographicCRS” property in your “projcrs” object might be due to the specific implementation or limitations of the “projcrs” class in your environment.

To retrieve the Earth radius of the Authalic sphere, you can consider an alternative approach. Instead of relying on the “GeographicCRS” property, you can directly access the projection parameters of the “projcrs” object. In your example, the Earth radius is specified as a parameter in the “WKT string”, so you can extract it from there.Here's an example of how you can extract the Earth radius from the WKT string using regular expressions in MATLAB:

 
wkt = "PROJCS[""MODIS Sinusoidal"",BASEGEOGCRS[""User"",DATUM[""World Geodetic Survey 1984"",SPHEROID[""Authalic_Spheroid"",6371007.181,0.0]],PRIMEM[""Greenwich"",0.0],UNIT[""Degree"",0.0174532925199433]],PROJECTION[""Sinusoidal""],PARAMETER[""False_Easting"",0.0],PARAMETER[""False_Northing"",0.0],PARAMETER[""Central_Meridian"",0.0],UNIT[""Meter"",1.0]]";
  
 
  
% Extract the Earth radius from the WKT string
  
radiusPattern = 'SPHEROID\[""Authalic_Spheroid"",(\d+\.\d+),0\.0\]';
  
radiusMatch = regexp(wkt, radiusPattern, 'tokens', 'once');
  
earthRadius = str2double(radiusMatch{1});
  
 
  
disp(earthRadius);

This code snippet uses regular expressions to extract the Earth radius value from the “WKT” string. The pattern “SPHEROID\[""Authalic_Spheroid"",(\d+\.\d+),0\.0\]” matches the part of the string that represents the Earth radius. The extracted value is then converted to a numeric type using “str2double()”.

By using this approach, you can obtain the Earth radius without relying on the “GeographicCRS” property of the “projcrs” object.

Here are some MATLAB documentations for you reference. Projected coordinate reference system object - MATLAB (mathworks.com) Convert strings to double precision values - MATLAB str2double (mathworks.com)

Hope it helps!