Issue on admin_level=4

Good morning all,

we tried to import argentina-latest.osb.pbf in a postgis db.
When we try to extract the data to create a layer of only boundaries with admin_level = ‘4’ with this query:

insert into osm_argentina.states(osm_id, name, uppername, geom)
SELECT planet_osm_polygon.osm_id,
planet_osm_polygon.name as name,
upper(planet_osm_polygon.name) AS uppername,
st_multi(planet_osm_polygon.way)::geometry(MultiPolygon, 3857) as way
FROM planet_osm_polygon
WHERE planet_osm_polygon.admin_level = ‘4’::text AND planet_osm_polygon.boundary = ‘administrative’::text;

We get the following rendering error, as you can see the boundary is not the geographic one, but the political.

Could it be that there are some error in the original .osm.pbf file?

Unfortunally this is not the only area of south America with this kind of issue, there are some other area which we cannot give an explanation

Is it a known issue? Is there a solution? Maybe download an older version of the area?

Thank you.

Hello too,

boundaries with an admin_level form the administrative level of a feature within a subdivision hierarchy in OpenStreetMap.

The geographical border, in other words presumably the land mass of a state, is only occasionally maintained in OpenStreetMap and does not appear to exist for Argentina.

Poland for example does have both:

As far as I can see there is no land_area relation at all in South America: overpass turbo

In the end and in other words: this is a feature and not a bug :slight_smile:

1 Like

Thank you for the explanation, so in other words there is no way to extract the land area?

Sorry but I can’t understand how openstreemap (or geofabrik for example) are able to render the landmass of the south america?

Another point that I do no agree is the fact that not all the political area is drawn, but only a portion in the Chubut area, and this is pretty strange (see purple line)

arg1

What do you mean by that? Can you provide a link or screenshot where the landmass can be seen?
By the way, admin_level is naturally an administrative boundary, so not sure what you are expecting.

Edit: There is the Cono Sur relation which maps the southern coastline of South America. You might have some success if you remove the areas that lie outside the Cono Sur.

1 Like

As @mcliquid wrote, the elements with

[code]
boundary=administrative
admin_level=4
[/code}

are borders that you would call “political” ones. The complete border of the state of Chubut would be the ways (not the nodes) of the relation Relation: ‪Chubut‬ (‪153548‬) | OpenStreetMap

The rendering of the coastline will be probably done with help of the ways that form the coastline

natural=coastline

like Relation: ‪Chubut‬ (‪153548‬) | OpenStreetMap These are elements of some relations like the one mentioned by @Jofban or

that you probably want to “remove” from the the admin_level=4-border of Chubut.

I’m sorry that I can’t help you with the needed SQL, but I hope I gave you some insight of the meaning of some of the elements in OpenStreetMap that you’re interested in.

Adding an independent land mass relation to OpenStreetMap is the wrong solution because it can be derived using geometric functions. You just have to intersect the land mass with the boundary polygon:

Download land mass polygons from osmdata.openstreetmap.de.

wget https://osmdata.openstreetmap.de/download/land-polygons-complete-3857.zip
unzip land-polygons-complete-3857.zip

Load the data into the PostgreSQL database using ogr2ogr (documentation of PostgreSQL driver):

PG_USE_COPY=YES ogr2ogr -f PostgreSQL PG:dbname=gis land-polygons-complete-4326/land_polygons.shp -nln land_polygons -gt 10000

Run intersection:

INSERT INTO osm_argentina.states (osm_id, name, uppername, geom)
SELECT
    p.osm_id AS osm_id,
    p.name AS name,
    UPPER(p.name) AS uppername,
    ST_Intersection(p.way, ST_Union(l.geom)) AS geom
  FROM planet_osm_polygon AS p
  JOIN land_polygons AS l ON (p.way && l.geom AND ST_Intersects(p.way, l.geom)
  WHERE p.admin_level = '4' AND p.boundary = 'administrative' AND p.osm_id < 0
  GROUP BY p.osm_id, p.admin_level, p.name;

For better performance:

  • Ensure that land_polygons has a spatial index.
  • Add a conditional index on the geometry column of planet_osm_polygon matching your WHERE condition: CREATE INDEX polygon_idx_admin4 ON planet_osm_polygon USING gist(way) WHERE admin_level = '4' AND boundary = 'administrative' AND osm_id < 0;
  • WHERE osm_id < 0 because boundary polygons are relations. Relations have negative IDs in the database if they are loaded with Osmpgsql’s old database layout as you do.
2 Likes

After importing the .osm.pbf file of south america and extracted the admin_level 4, the result of the boundaries is the following

As you can see the real coastline (geographical line) is replaced by what I think is the political line, but this happens only in Chubut and some other area

Ok thank you, I’ll try immediately

I can’t tell you why the mappers decided to map differently between different states like

and Chubut but you should never assume that an administrative border ends at a coastline.

@Nakaner wrote how to subtract the maritime parts from the administrative parts.

Please be aware that in OpenStreetMap boundary=administrative is not the same as boundary=political :slight_smile:

I have similar problem also in my forestpark layer

drop table if exists forestpark;
create table forestpark(
  id serial not null primary key,
  osm_id integer,
  name text,
  geom geometry(multipolygon, 3857)
);
create index gix_forestpark on forestpark using gist(geom);
delete from forestpark;
insert into forestpark(osm_id, name, geom) 
  SELECT planet_osm_polygon.osm_id, 
    planet_osm_polygon.name, 
    st_multi(planet_osm_polygon.way)::geometry(MultiPolygon, 3857) as way
  FROM planet_osm_polygon
  WHERE (planet_osm_polygon.natural = 'wood') OR (planet_osm_polygon.landuse = ANY (ARRAY['greenery'::text, 'green'::text, 'garden'::text, 'Reserve_forest'::text, 'forest'::text, 'orchard'::text, 'park'::text, 'plant_nursery'::text, 'grass'::text, 'greenfield'::text, 'meadow'::text, 'recreation_ground'::text, 'village_green'::text, 'vineyard'::text])) OR (planet_osm_polygon.leisure = ANY (ARRAY['nature_reserve'::text, 'park'::text, 'dog_park'::text, 'garden'::text, 'golf_course'::text, 'horse_riding'::text, 'recreation_ground'::text, 'stadium'::text]));
--delete from forestpark where not st_intersects(st_centroid(geom), (select geom from country limit 1));

Could you please fix the query?

Have you tried looking for yourself at what is in the OSM data in the areas causing you problems?

For example, it’s fairly obvious that the large area south of South America is this marine nature reserve:

If you look at the tags on that object, and compare with the tags in your query, you should see why it is included in your results.

“Fixing” the query would then depend on what exactly you want in your “forestpark”. The query seems to include a long list of tags, some of which represent quite different objects. If you want all nature reserves (among many other things) there may be nothing to fix.

3 Likes