0

In code one below, there a query that finds the points of intersection in terms of longitude and latitude. And in code two, it is to show specific info like type, properties and geometry, etc.

What I want to achieve is to have the main query in code two contains information about the point of intersection "longitude and latitude" and the area as well. In other words, given the geom in code two, I would like to integrate code one into code two to have information about the points of intersection and area.

How can code one be integrated into code two?

code1_to find the coordinates of intersection:

query ="""SELECT ST_X(ST_Transform(point,4326)) as lon, ST_Y(ST_Transform(point,4326)) as lat, ST_AsText(ST_Transform(point,4326)),ST_Area(
                ST_Intersection(
                    ST_SetSRID(
                        ST_MakeEnvelope(ST_X(point),ST_Y(point),ST_X(point)+{width}, ST_Y(point)+{height}),25832),
                            ST_Transform(
                                    ST_SetSRID(ST_GeomFromGeoJSON(
                                        '{geometry}'),4326)
                                        ,25832)))
    FROM {table} 
    WHERE 
    st_intersects(
        ST_Transform(
            ST_SetSRID(ST_GeomFromGeoJSON(
               '{geometry}'),4326)
               ,25832),
                st_setsrid(ST_MakeEnvelope(st_x(point),st_y(point),st_x(point)+{width},st_y(point)+{height}),25832))""".format(table=config['PostgreDB']['table_name_test'], width=config['Grid']['cell_width'], height=config['Grid']['cell_height'],geometry=geometry)        
                

code2:

query = """  WITH data AS (
        SELECT '{featuresCollection}'::json AS featuresCollection
        )
        SELECT gid,geom,type::text,properties::text,
        array_to_string(array_agg(x_4326||' '||y_4326 ORDER BY gid),',') AS g4326,
        array_to_string(array_agg(x_25832||' '||y_25832 ORDER BY gid),',') AS g25832             
        FROM (
        SELECT
        ROW_NUMBER() OVER () AS gid,
        ST_AsText(ST_GeomFromGeoJSON(feature->>'geometry')) AS geom,
        feature->>'type' AS type,
        feature->>'properties' AS properties,
        ST_X((ST_DumpPoints((ST_GeomFromGeoJSON(feature->>'geometry')))).geom) x_4326,       
        ST_Y((ST_DumpPoints((ST_GeomFromGeoJSON(feature->>'geometry')))).geom) y_4326,  
        ST_X((ST_DumpPoints((ST_Transform(ST_SetSRID(ST_GeomFromGeoJSON(feature->>'geometry'),4326),25832)))).geom) x_25832,       
        ST_Y((ST_DumpPoints((ST_Transform(ST_SetSRID(ST_GeomFromGeoJSON(feature->>'geometry'),4326),25832)))).geom) y_25832       

        FROM (SELECT json_array_elements(featuresCollection->'features') AS feature FROM data) AS f) j
        GROUP BY gid,type::text,properties::text,geom
        ORDER BY gid;""".format(featuresCollection=featuresCollection)

Sample data:

[(3338490, 5668960, Decimal('1.02'), Decimal('52.08'), '0101000020E864000077D23C26C5A81441A9BAEC5A4F9E5541'), (3338490, 5668950, Decimal('0.77'), Decimal('52.13'), '0101000020E864000047A52726C5A81441D4552EDB4C9E5541'), (3338490, 5668940, Decimal('0.36'), Decimal('52.19'), '0101000020E864000005781226C5A8144109F16F5B4A9E5541')]

Image of some data in the table:

enter image description here

Data in the featureCollection:

DB Fiddle

9
  • could you also provide some sample data and show the exact expected result? Commented May 26, 2021 at 12:39
  • @JimJones i added a link to the data.it is the same data as in the previous one.please have a look. the data in the first select statement:SELECT '{"type": "FeatureCollection", "features": [{"type": "Feature", "properti Commented May 26, 2021 at 12:47
  • Can you also add the exact expected result? I'm confused about code 1. Where does the variable point come from? Where is the GeoJSON for {geometry}? Commented May 26, 2021 at 12:57
  • @JimJones for point it is a column in table "s_30m_test" for '{geometry}' it must be substituted with geom from "feature->>'geometry'" Commented May 26, 2021 at 13:06
  • Can you add some data sample from table s_30m_test and the exact expected result based on the given json document? Commented May 26, 2021 at 13:07

1 Answer 1

0

Just place the query in code 2 in the FROM clause and either join it with code 1 or just match them in the WHERE clause, e.g.

    query = """ 
        SELECT j.*,
            ST_X(ST_Transform(point,4326)) As lonOfIntersection, 
            ST_Y(ST_Transform(point,4326)) AS latOfIntersection, 
            ST_AsText(ST_Transform(point,4326)) pointOfIntersectionEPSG4326,
            ST_AsText(ST_Transform(point,25832)) pointOfIntersectionEPSG25832,
            ST_Area(
                ST_Intersection(
                ST_SetSRID(
                    ST_MakeEnvelope(
                    ST_X(point),
                    ST_Y(point),
                    ST_X(point)+{width}, 
                    ST_Y(point)+{height}),
                    25832),j.geometry
                )
            ) As areaOfCoverage
        FROM {table}
        JOIN (
            WITH data AS (
            SELECT '{featuresCollection}'::json AS featuresCollection
            )
            SELECT DISTINCT
            geometryID,geomType,geomProperties,
            array_to_string(array_agg(x_4326||' '||y_4326 ORDER BY geometryID),',') AS polygonsCoordinatesInEPSG4326,
            array_to_string(array_agg(x_25832||' '||y_25832 ORDER BY geometryID),',') AS polygonsCoordinatesInEPSG258,
            geometry
            FROM (
            SELECT 
                ROW_NUMBER() OVER () AS geometryID,
                feature->>'type' AS geomType,
                feature->>'properties' AS geomProperties,
                ST_X((ST_DumpPoints((ST_GeomFromGeoJSON(feature->>'geometry')))).geom) AS x_4326,       
                ST_Y((ST_DumpPoints((ST_GeomFromGeoJSON(feature->>'geometry')))).geom) AS y_4326,   
                ST_X((ST_DumpPoints((ST_Transform(ST_SetSRID(ST_GeomFromGeoJSON(feature->>'geometry'),4326),25832)))).geom) AS x_25832,       
                ST_Y((ST_DumpPoints((ST_Transform(ST_SetSRID(ST_GeomFromGeoJSON(feature->>'geometry'),4326),25832)))).geom) AS y_25832,
                ST_Transform(ST_SetSRID(ST_GeomFromGeoJSON(feature->>'geometry'),4326),25832) AS geometry
 
            FROM (SELECT json_array_elements(featuresCollection->'features') AS feature 
                    FROM data) AS f) j
            GROUP BY geometryID,geometry,geomType,geomProperties) j ON   
            ST_Intersects(j.geometry,
                ST_SetSRID(
                    ST_MakeEnvelope(
                    ST_X(point),
                    ST_Y(point),
                    ST_X(point)+{width},
                    ST_Y(point)+{height}),25832));
""".format(table=config['PostgreDB']['table_name_test'], width=config['Grid']['cell_width'], height=config['Grid']['cell_height'],featuresCollection=featuresCollection)
Sign up to request clarification or add additional context in comments.

3 Comments

Hi, would you like to answer this question:stackoverflow.com/questions/67738265/…
Morning, i would test and acknowledge regarding your answer on Monday...have a nice weekend
hi, would you please answer this question:stackoverflow.com/questions/68558791/…

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.