[PROJ] Affine transformation using DERIVEDPROJCRS
Even Rouault
even.rouault at spatialys.com
Fri Jul 8 13:29:17 PDT 2022
Bart,
> ** RE-SENT AS PLAIN TEXT IN VIEW OF EMAIL SIZE **
>
> Hi Even,
>
> I worked things through, and have included the following three routines in test_crs.cpp, that compile and test okay.
> I have a few questions / remarks though:
> 1. Why is the EPSG definition ID of an affine transform not being exported: ID["EPSG",9624]],
You need to add it to the PropertyMap with
.set(metadata::Identifier::CODESPACE_KEY, metadata::Identifier::EPSG)
.set(metadata::Identifier::CODE_KEY, 9624)
> 2. Same for the ID’s of the Affine transformation parameters: 8623, 8624, 8625, 8639, 8640 and 8641.
Similarly as above on the property map of each parameter.
> 3. The USAGE description isn’t being included in the generated WKT string and is also missing in the derived transform
> USAGE[
> SCOPE["Engineering survey, topographic mapping."],
> AREA["Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone."],
> BBOX[50.75,3.2,53.7,7.22]]]
> This is really a pity. Is there an ‘easy’ way to fix this ?
There's no reason that a DerivedProjectedCRS should have the same usage
as its base one.
You can re-use the one of the base CRS (the projected CRS) with
something along
propertyMap.set(ObjectUsage::OBJECT_DOMAIN_KEY, baseCRS->domains()[0]);
If you want to build a custom objectDomain object, the best example
would be to look at the JSONParser::buildObjectDomain() implementation
in src/iso19111/operation/oputils.cpp
where propertyMap is the first argument of
DerivedProjectedCRS::create(), on which you'll also call
.set(IdentifiedObject::NAME_KEY, "Abtswoudsebrug")
> 4. Nevertheless, maybe the example below is more ‘meaningful’ than the current one, that is included in the test_crs.file.
Those unit tests are just that, unit tests, so they indeed tend to be
minimalist.
Perhaps a good contribution would be an Example section under
https://proj.org/development/index.html
Even
>
> Kind regards,
>
> Bart
>
>
> ProjectedCRSNNPtr createProjected_v2()
> {
> auto dbContext = DatabaseContext::create();
> auto factory = AuthorityFactory::create(dbContext, "EPSG");
> auto crs = factory->createCoordinateReferenceSystem("28992");
> auto proj_crs = nn_dynamic_pointer_cast<ProjectedCRS>(crs);
> EXPECT_TRUE(proj_crs != nullptr);
>
> return NN_CHECK_ASSERT(proj_crs);
> }
>
> static DerivedProjectedCRSNNPtr createDerivedProjectedCRS_v2() {
>
> auto derivingConversion = Conversion::create
> (
> PropertyMap().set(IdentifiedObject::NAME_KEY, "Affine"),
> PropertyMap().set(IdentifiedObject::NAME_KEY, "Affine parametric transformation"),
> std::vector<OperationParameterNNPtr>
> {
> OperationParameter::create(PropertyMap().set(
> IdentifiedObject::NAME_KEY, "A0")),
> OperationParameter::create(PropertyMap().set(
> IdentifiedObject::NAME_KEY, "A1")),
> OperationParameter::create(PropertyMap().set(
> IdentifiedObject::NAME_KEY, "A2")),
> OperationParameter::create(PropertyMap().set(
> IdentifiedObject::NAME_KEY, "B0")),
> OperationParameter::create(PropertyMap().set(
> IdentifiedObject::NAME_KEY, "B1")),
> OperationParameter::create(PropertyMap().set(
> IdentifiedObject::NAME_KEY, "B2")),
> },
> std::vector<ParameterValueNNPtr>
> {
> ParameterValue::create(Length(-257762.015)),
> ParameterValue::create(Scale ( 0.914845326707141)),
> ParameterValue::create(Scale ( 0.403804443019272)),
> ParameterValue::create(Length(-374179.972)),
> ParameterValue::create(Scale (-0.403804443019272)),
> ParameterValue::create(Scale ( 0.914845326707141)),
> }
> );
>
> // create a DerivedProjectedCRS, but ...
> // first we need to get the BaseCRS from an EPSG code.
>
> return DerivedProjectedCRS::create(
> PropertyMap().set(IdentifiedObject::NAME_KEY, "Abtswoudsebrug"),
> createProjected_v2(), derivingConversion,
> CartesianCS::createEastingNorthing(UnitOfMeasure::METRE));
> }
>
> TEST(crs, derivedProjectedCRS_WKT2_2019_v2) {
>
> auto expected_v2 =
> "DERIVEDPROJCRS[\"Abtswoudsebrug\",\n"
> " BASEPROJCRS[\"Amersfoort / RD New\",\n"
> " BASEGEOGCRS[\"Amersfoort\",\n"
> " DATUM[\"Amersfoort\",\n"
> " ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,\n"
> " LENGTHUNIT[\"metre\",1]]],\n"
> " PRIMEM[\"Greenwich\",0,\n"
> " ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
> " CONVERSION[\"RD New\",\n"
> " METHOD[\"Oblique Stereographic\",\n"
> " ID[\"EPSG\",9809]],\n"
> " PARAMETER[\"Latitude of natural origin\",52.1561605555556,\n"
> " ANGLEUNIT[\"degree\",0.0174532925199433],\n"
> " ID[\"EPSG\",8801]],\n"
> " PARAMETER[\"Longitude of natural origin\",5.38763888888889,\n"
> " ANGLEUNIT[\"degree\",0.0174532925199433],\n"
> " ID[\"EPSG\",8802]],\n"
> " PARAMETER[\"Scale factor at natural origin\",0.9999079,\n"
> " SCALEUNIT[\"unity\",1],\n"
> " ID[\"EPSG\",8805]],\n"
> " PARAMETER[\"False easting\",155000,\n"
> " LENGTHUNIT[\"metre\",1],\n"
> " ID[\"EPSG\",8806]],\n"
> " PARAMETER[\"False northing\",463000,\n"
> " LENGTHUNIT[\"metre\",1],\n"
> " ID[\"EPSG\",8807]]]],\n"
> " DERIVINGCONVERSION[\"Affine\",\n"
> " METHOD[\"Affine parametric transformation\"],\n"
> " PARAMETER[\"A0\",-257762.015,\n"
> " LENGTHUNIT[\"metre\",1,\n"
> " ID[\"EPSG\",9001]]],\n"
> " PARAMETER[\"A1\",0.914845326707141,\n"
> " SCALEUNIT[\"unity\",1,\n"
> " ID[\"EPSG\",9201]]],\n"
> " PARAMETER[\"A2\",0.403804443019272,\n"
> " SCALEUNIT[\"unity\",1,\n"
> " ID[\"EPSG\",9201]]],\n"
> " PARAMETER[\"B0\",-374179.972,\n"
> " LENGTHUNIT[\"metre\",1,\n"
> " ID[\"EPSG\",9001]]],\n"
> " PARAMETER[\"B1\",-0.403804443019272,\n"
> " SCALEUNIT[\"unity\",1,\n"
> " ID[\"EPSG\",9201]]],\n"
> " PARAMETER[\"B2\",0.914845326707141,\n"
> " SCALEUNIT[\"unity\",1,\n"
> " ID[\"EPSG\",9201]]]],\n"
> " CS[Cartesian,2],\n"
> " AXIS[\"(E)\",east,\n"
> " ORDER[1],\n"
> " LENGTHUNIT[\"metre\",1,\n"
> " ID[\"EPSG\",9001]]],\n"
> " AXIS[\"(N)\",north,\n"
> " ORDER[2],\n"
> " LENGTHUNIT[\"metre\",1,\n"
> " ID[\"EPSG\",9001]]]]";
> // " USAGE[\",\n"
> // " SCOPE[\"Engineering survey, topographic mapping.\"],\n"
> // " AREA[\"Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone.\"],\n"
> // " BBOX[50.75, 3.2, 53.7, 7.22]]]";
>
> auto crs = createDerivedProjectedCRS_v2();
> EXPECT_EQ(
> crs->exportToWKT(
> WKTFormatter::create(WKTFormatter::Convention::WKT2_2019).get()), expected_v2);
>
> EXPECT_TRUE(crs->isEquivalentTo(crs.get()));
> EXPECT_TRUE(crs->shallowClone()->isEquivalentTo(crs.get()));
> EXPECT_FALSE(crs->isEquivalentTo(createUnrelatedObject().get()));
>
> auto geodCRS = crs->extractGeodeticCRS();
> EXPECT_TRUE(geodCRS != nullptr);
> }
>
>
>
--
http://www.spatialys.com
My software is free, but my time generally not.
More information about the PROJ
mailing list