texts: How to speed up raster/map algebra in PostGIS?

texts:
I am testing map algebra with PostGIS (2.1 with PostgreSQL 9.3) on two images, and found it to be rather slow. For example, dividing two images of 7200 * 3600 pixels took 110 ~ 120s on a desktop with i7 3.4 GHz CPU; the same operation with the same data took ArcGIS (10.1) 15 ~ 16s.

The SQL statement tested was based on this PostGIS example as follows:

SELECT ST_MapAlgebra(t1.rast,t2.rast
        ,'([rast1]+0.001)/([rast2]+0.001)') AS rast
FROM img1 t1, img2 t2

The two rasters were imported using the following options:

 raster2pgsql -d -M -C -I -s 4326 ...

I did not use filesystem storage or tiling.

My question is:

What's the best way (if possible) to improve the above SQL for map algebra?


texts: Will ArcGIS for Desktop run on Windows 10?

texts:
Will ArcGIS 10 run on the new windows 10 technical preview?

I am setting up an interface on my Mac and there is the free option of downloading the Windows 10 technical preview, some sort of beta version, I'm assuming.

Just wanted to see if 10.1 will run ok on Windows 10 (ignoring the Mac aspect), as this is all I am using my windows interface for?


texts:
I have developed an AddIn for ArcGIS 10.4 using Visual Studio 2015. My AddIn has several "buttons" defined in the config.esriaddinx XML file.

It is possible to completely hide a button when ArcMap starts? I know that I can override the button OnUpdate() event to disable the command, but I want to make it invisible to users.

My AddIn does have an Extension and I'm hoping I can put some code in the extension OnStartUp() event to either hide or remove a button.


texts: Extending batch of lines to another object using ArcGIS Desktop?

texts:
Is there any simple tool which allows to extend all lines in .shp to another .shp?

To clarify, I do have river.shp and sections.shp on the river. I'd like to produce buffer from river, extend lines in section.shp to the buffer or clip lines if neccesarly (may happend in some cases).

I thought it would be a simple operation, however can't find any tool which would extend my lines. Do I miss something? Any suggestions?


texts: Rasterizing polygons results in gaps that don't appear in the original polygons

texts:
I'm getting gaps when rasterizing polygons and was hoping to get an idea on the best solution.
The first picture is the polygons and there doesn't appear to be any visible gaps on the edges of the polygons even when zoomed in as close as possible on QGIS.
The second picture is the rasterized versions and you can see a gap between the two rasters. This is in UTM18N around NY and NJ.


texts: Removing Black Edge in Aerial Imagery using ArcGIS Pro

texts:
How do I remove the black edge in image below on a raster aerial imagery?
I am using ArcGIS Pro.


texts:
Here is the "big" problem: my Digital elevation model is divided in a lot of pieces (single files .xyz).

In order to build a mosaic of them I need to import them all. With SAGA I can import them just one by one (Import grid from xyz). Is there an option or an easier way to import them all together?

P.S.: I'm using SAGA 2.1.2


texts: How to export/share my personalized styles on QGIS?

texts: I created a lot of personalized styles for QGIS and my co-workers are asking me to share the styles with them. How can I export my personalized styles of lines, areas and points?


texts: Issues creating slope raster layer from a point layer

texts:
I'm trying to create a slope raster layer using QGIS 2.0 on Windows XP.

I started with a points layer with the Z elevation attribute.

I then used the interpolation plugin to create a DEM.

But, when I used the Raster Terrain Analysis plugin, it gave me a strange output. What can I be doing wrong?


texts: Append tool snippet not working in Pyscripter?

texts:
I need to append two shapefiles together. The append tool works within the ArcMap toolboxes AND Python window, but it does not work when copying over the python snippet to PyScripter. It'll completely run through but won't produce values within the fields (all values will be empty or null).

I'm assuming it has something to do with the FieldMapping() objects?

import arcpy

arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:\Users\Documents\ArcGIS\Default.gdb"

tc = r"testing_append"
changesshp = r"C:\Users\gfrederick\Documents\Scripts\changes.shp"

lyr1 = arcpy.MakeFeatureLayer_management(tc, "tc_lyr")
lyr2 = arcpy.MakeFeatureLayer_management(changesshp, "change.lyr")

arcpy.Append_management(lyr2, lyr1, schema_type="NO_TEST", field_mapping="""RECORDNUM "RECORDNUM" true false false 19 Double 10 18 ,First,#,changes_new,RECORDNUM,-1,-1;TAKENBY "TAKENBY" true false false 50 Text 0 0 ,First,#,changes_new,TAKENBY,-1,-1;SETDATE "SETDATE" true true false 8 Date 0 0 ,First,#,changes_new,SETDATE,-1,-1;SETYEAR "SETYEAR" true false false 4 Text 0 0 ,First,#;MCD "MCD" true false false 50 Text 0 0 ,First,#,changes_new,MCD,-1,-1;RDPREFIX "RDPREFIX" true false false 50 Text 0 0 ,First,#,changes_new,RDPREFIX,-1,-1;ROUTE "ROUTE" true false false 19 Double 10 18 ,First,#,changes_new,ROUTE,-1,-1;ROAD "ROAD" true false false 50 Text 0 0 ,First,#,changes_new,ROAD,-1,-1;CNTDIR "CNTDIR" true false false 50 Text 0 0 ,First,#,changes_new,CNTDIR,-1,-1;FROMLMT "FROMLMT" true false false 50 Text 0 0 ,First,#,changes_new,FROMLMT,-1,-1;TOLMT "TOLMT" true false false 50 Text 0 0 ,First,#,changes_new,TOLMT,-1,-1;TYPE "TYPE" true false false 20 Text 0 0 ,First,#,changes_new,TYPE,-1,-1;RDSUFFIX "RDSUFFIX" true false false 15 Text 0 0 ,First,#,changes_new,RDSUFFIX,-1,-1;LATITUDE "LATITUDE" true false false 19 Double 8 18 ,First,#,changes_new,LATITUDE,-1,-1;LONGITUDE "LONGITUDE" true false false 19 Double 8 18 ,First,#,changes_new,LONGITUDE,-1,-1;FACTOR "FACTOR" true false false 19 Double 8 18 ,First,#,changes_new,FACTOR,-1,-1;AXLE "AXLE" true false false 19 Double 8 18 ,First,#,changes_new,AXLE,-1,-1;AADT "AADT" true false false 9 Long 0 9 ,First,#,changes_new,AADT,-1,-1;UPDATED "UPDATED" true true false 8 Date 0 0 ,First,#;CO_NAME "CO_NAME" true false false 35 Text 0 0 ,First,#;MUN_NAME "MUN_NAME" true false false 40 Text 0 0 ,First,#""", subtype="")


texts: Split a polygon by another polygon and remove the overlapping area

texts: Using DotSpatial I want to split a polygon by another overlapping polygon. The overlapped area can be clipped using the Intersection method but I did not find any way to split the original source polygon and remove the overlapping area.


texts: What would you call the polygon corresponding to a coastline theme? Land mass? It corresponds to continents and islands, but I'm not sure about the best word. The opposite of the sea.


texts: unable to render/view a postgis raster in qgis 2.14.1

texts:
I have successfully added a raster to postgis using the following:

raster2pgsql -s 4326 -I -C -M 00N_030E.tif -F -t 250x250 -a test.00N_030E | psql -U postgres -d raster -h localhost -p 5432

I tried loading the raster into QGIS 2.14.1 to view it. However, it is listed as a table without a geometry.

This is from the QGIS message log:
2016-05-09T12:35:34    1   Field rast ignored, because of unsupported type raster

I have tried to enable all raster drivers using  the following:

ALTER DATABASE raster SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';

and

ALTER SYSTEM SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';

but did not seems to work.

Any suggestions as to what the issue maybe and how to rectify it?


texts: Set display name of vector layer in processing script

texts:
I have a processing script  (QGIS 3.12) that extracts values from an input table using the native:extractbyattribute algorithm. This produces a temporary output table with the display name 'Extracted (attribute) and I would like to rename this to something more meaningful within the script.

I can get the unique name of the layer with:

extracted = processing.run('native:extractbyattribute', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
layer = QgsProcessingUtils.mapLayerFromString(format(extracted['OUTPUT']), context )

but I can't see how to set the display name programmatically within the script. The layer is of type QgsVectorLayer. I can get its attributes with layer.dataProvider(), but this does not include a name attribute.

There are several ways to iterate through the layers from the console, but no information on how to do this through a script.


texts: Is there an automatic function in QGIS to package all files in a project for transferability?

texts:
In my project, I use a lot of geodata from different directories.
Now I want to give the QGIS project to a colleague (eg on a CD).

Is it possible to copy the shapes from different directories automatically to one directory?


texts: Install gdal2tiles inside conda environment

texts:
I need to use gdal2tiles with Anaconda 3 on Ubuntu 18.04. I've created the environment using:

  conda create --name gisenv python=3.6

I've installed GDAL with:

  conda install gdal

After this I've created a simple python script import gdal2tiles but I see this error message:

  ModuleNotFoundError: No module named 'gdal2tiles'

I thought that using conda install gdal2tiles would solve the problem but I see this message:

  Collecting package metadata (current_repodata.json): done Solving
  environment: failed with initial frozen solve. Retrying with flexible
  solve. Collecting package metadata (repodata.json): done Solving
  environment: failed with initial frozen solve. Retrying with flexible
  solve.

  PackagesNotFoundError: The following packages are not available from
  current channels:

  gdal2tiles

  Current channels:

  https://repo.anaconda.com/pkgs/main/linux-64
  https://repo.anaconda.com/pkgs/main/noarch
  https://repo.anaconda.com/pkgs/r/linux-64
  https://repo.anaconda.com/pkgs/r/noarch

  To search for alternate channels that may provide the conda package
  you're looking for, navigate to

https://anaconda.org

  and use the search bar at the top of the page.

How I can solve?


texts:
I want to restrict the OpenLayers 3 heatmap-earthquakes examle to show only Europe.

I have:

var raster = new ol.layer.Tile({
  source: new ol.source.Stamen({
    layer: 'toner-background'
  })
});

var map = new ol.Map({

  layers: [raster, vectorLayer],
  target: 'map',
    view: new ol.View({
        projection: 'EPSG:3857',
        center: [0, 0],
        maxZoom: 9,
        minZoom: 0,
        zoom: 4
  })
});

Let`s say that the points:

var point1 = new ol.geom.Point(ol.proj.transform([-20.000000, 28.000000], 'EPSG:4326', 'EPSG:3857'));
var point2 = new ol.geom.Point(ol.proj.transform([50.000000, 73.000000], 'EPSG:4326', 'EPSG:3857'));

are the bounds of Europe, how can I restrict the map view only between the two point?


texts:
I want to create mbtiles with the GDAL Python API. How can I integrate format specific options like "TILE_FORMAT = 'jpeg', QUALITY='90'" for mbtiles?
from osgeo import gdal

path = r'C:\Users\go\my.vrt'
pathout = r'C:\Users\go\my.mbtiles'
src_ds = gdal.Open(path)
ds = gdal.Translate(pathout , src_ds, format = 'mbtiles')
ds = No


texts: Splitting line into equal area segments and assigning sequential value from start to end

texts:
I have a line file that I would like to divide into 1-km segments and assign each segment a sequential value from the beginning of the line (1) to the end (total number of 1-km sections. I found someway of doing this in another program (FME), but I am wondering if there is a way of doing this in ArcGIS or possibly R.

In Arc, I have been working with the COGO proportion tool. This does a nice job of splitting the line into segments, but I can't figure out how to assign these new segments the sequential values from start to end.


texts: call child div function from parent (not using window.opener)

texts:
This is a a Javascript/JQuery question that is directly related to the ArcGIS Javascript API...

I have a map div created using a 'data' tag in HTML, which becomes visible on a button click.

(I'll call the map the 'child' and main page the 'parent' although I am not sure that is technically correct in this case?)

I would like to force a resize of the map when a 'fullscreen' button is clicked in the parent page.  Generally the map resizes on its own but sometimes it does not...

How can I do this? I did not create a new window with 'window.open' ...

My HTML for the child window is:

<div style='width:100%;  height:100%'><object type="text/html" style='width:100%;  height:100%' data='mapPage.html'>

..and the JS for the fullscreen button :

    mapFullscreen.click(function () {
    window.child.map.refresh(); // this doesnt seem to do anything

    $(mapContainer).animate({
        "height": $(window).outerHeight(),
        "width": $(window).outerWidth(),
       "top": 0,
       "left": 0

   }, 300, function(){
    $(".map-fullscreen").fadeOut();
  });
});

Thanks!


texts: Raster - attribute table, fields - confused

texts:
I am a bit confused, as during my studies I learnt that rasters have only one cell value and for some it might be true, as the attribute table is greyed out in ArcGIS. But for others, multiple fields are provided, so the assumption seems to be wrong.

Can anyone explain me the difference of rasters that can contain one or multiple fields?

In ArcGIS an attribute table can be created for raster, but only for integer raster, not for float raster (see the "Build Raster Attribute Table" function).


texts: How to add GeoJSON layer with GraphHopper?

texts:
I have a GeoJSON file containing some POI that I would like to be able to display within a separate layer with GraphHopper.
But after several tries and search over internet, i can't manage to get a way to do it.

The idea is to display additional GeoJSON layers over the built-in GraphHopper ones.

Here is the code, https://graphhopper.com/#community

and what i am trying to do is to get an additional layer to display POI information from a GeoJSON file.

Here is a sample of the information available in the GeoJSON file :

{"type": "Feature",
"properties": {
    "fee": "no",
    "bicycle_parking": "anchors",
    "ref": "PVNAN23",
    "address": "Rue Gabriel Goudy 44200 Nantes",
    "name": "Pirmil P+R",
    "capacity": "24",
    "park_ride": "yes",
    "amenity": "bicycle_parking",
    "covered": "yes"
},
"geometry": {"type": "Point", "coordinates": [-1.5406709, 47.1960031]}},
{"type": "Feature",
"properties": {
    "bicycle_parking": "stands",
    "addr:postcode": "44000",
    "addr:country": "FR",
    "name": "Madeleine",
    "capacity": "6",
    "amenity": "bicycle_parking",
    "addr:street": "chaussée de la Madeleine",
    "note": "vérifié",
    "addr:city": "Nantes",
    "covered": "no",
    "addr:housenumber": "35"
},
"geometry": {"type": "Point", "coordinates": [-1.55076671448, 47.21000114109]}}

]}

The whole file has been validated in JSON validator and found no error.

I tried what is explained in this thread but nothing is working.
How to load external GeoJSON file into Leaflet map


texts: How to use map composer in a stand-alone script?

texts:
I'm trying to follow the map rendering section from the pyqgis cookbook, but I'd like to test this as a standalone application. I can do the first part, using simple rendering, but I'm a bit stuck doing the second example using the map composer as a standalone script.
Here is a standalone example for the bit I can do:
from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtXml import *

QgsApplication.setPrefixPath("/usr/", True)
QgsApplication.initQgis()

fh = open("eg.csv","w")
fh.write("""
x,y,name
153.0278, -27.4679, Brisbane
144.2500, -23.4500, Longreach
145.7753, -16.9256, Cairns
""")
fh.close()

uri = "eg.csv?delimiter=%s&xField=%s&yField=%s" % (",", "x", "y")
layer = QgsVectorLayer(uri, "eglayer", "delimitedtext")
QgsMapLayerRegistry.instance().addMapLayer(layer)
img = QImage(QSize(800,600), QImage.Format_ARGB32_Premultiplied)
color = QColor(255,255,255)
img.fill(color.rgb())
p = QPainter()
p.begin(img)
render = QgsMapRenderer()
lst = [ layer.getLayerID() ]  # add ID of every layer
render.setLayerSet(lst)
rect = QgsRectangle(render.fullExtent())
rect.scale(1.1)
render.setExtent(rect)
render.setOutputSize(img.size(), img.logicalDpiX())
render.render(p)
p.end()
img.save("render.png","png")

What I'd really like to do is the same, but use QgsComposition, and save as for example pdf. The cookbook says:

When using composer in a standalone application, you can create your own map renderer instance the same way as shown in the section above and pass it to the composition.

This bit I couldn't do, all my attempts either get an empty map, or a segfault.
I'm running linux mint 13, using qgis 1.8.0. It would be great if someone could show me how to modify the simple example to one that uses the composer.


texts: Weighting amount of overlapping polygons in postgis

texts:
I have two tables that I join in my existing workflow. One of these includes flood events (was a shapefile that has been imported to PostgreSQL), the other one is a grid that holds the amount of precipitation for each day and each grid. It is a grid spawned across the world, each pattern one degree high and one degree wide (64800 in total).
In my workflow, I combine the flood event with the summed up precipitation in the affected area.
And this is how I do it: At first I calculate the centroid of each pattern of the grid (st_centroid).

Then I check if the centroid is within an event of the flood-polygon, using the st_within function:

 CREATE TABLE within
    SELECT  cell_id, (id of the grid)
            geometry AS cell_geom,
            gid, (id of the flood event)
            started_at,
            finished_at,
            geom AS dfo_geom
    FROM    precipitation_grid,
            flood_event_lib
    WHERE
          ST_Within(precipitation_grid.centroid, flood_event_lib.geom);

CREATE INDEX within_gix ON within USING GIST (dfo_geom);

In the next step, I combine the two libraries, based on the cell_id that I used in the first step.

CREATE TABLE libjoin AS
    SELECT cell_id,
           gid AS dfo_id,
           started_at,
           finished_at,
           date,
           precipitation,
           flood_lib_geom,
           cell_geom
    FROM within, precipitation_lib
    WHERE within.cell_id = precipitation_lib.cell AND
(precipitation_lib.date >= within.started_at AND precipitation_lib.date <= within.finished_at)
    AND precipitation > 0

ORDER BY gid, cell_id, date;

In the final two steps, I calculate first the average of the daily precipitation

 CREATE TABLE libjoin_avg AS
   SELECT *,
   AVG(precipitation) OVER (PARTITION BY date ORDER BY dfo_id)
   FROM libjoin

And then I sum it up by adding the averages of each day within the flood event

 CREATE TABLE libjoin_sum AS
   SELECT *,
   SUM(avg) OVER (PARTITION BY started_at ORDER BY dfo_id)
   FROM libjoin_avg

This workflow is working but I would like to optimise it. My grid is not very precise and, as I am using the centroid of a pattern, I miss pretty many flood events because the centroid of a pattern is not always within the flood event. Another problem is, that once the centroid is within the flood event, the pattern is used for the calculation, no matter how much of its area is actually covered. You can see that in the following image, especially on the upper right side. These patterns are used for calculation because the centroid is within the flood event geometry.

This image shows my different layers, imported in QGIS. The green background is the actual grid and the highlighted blue ones are those with the centroid within the flood event. So they are used for the calculation of the summarised precipitation of the flood event shown in light blue on top.

Do you have any ideas of how I could optimise my workflow? I would like to weight the patterns by the amount of area that is covered by the flood event polygon and use that weighted amount for the calculation of the overall precipitation for each flood event but I don't know how to do that.

I am open for other solutions, not based on the centroid of the pattern, as well.


texts: How to select and delete past freehand features in ArcMap?

texts:
I am editing an older map in ArcMap 10.1 and want to edit a freehand drawing someone made. I cannot "select" it or find the shapefile it was made under...

Any advice?

I know how to make my own freehand drawing, but I also need to delete/edit the original one.


texts: How can I resample a image or imagecollection to a higher resolution in Google Earth Engine Python API

texts:
I try to use resample and reproject function in GEE to downscale the MODIS 250 products into 240m in order to correspond the 8*8 Landsat 30m pixels for furture purpose.

I used a Landsat image for example. Change image resolution in GEE is implemented by ee.image.reproject() function. Followed by the resolution given by Kersten and Tyler Erickson, i re-wrote the JS code into Python as follows:

import ee
ee.Initialize()

Reference a Landsat scene.

image_30m = ee.Image('LANDSAT/LC08/C01/T1/LC08_044034_20170614')

get the image projection by querying the projection of band2

band2 = image_30m.select('B2')
proj = band2.projection().getInfo()
crs = proj['crs']

resample and reproject the image to 10m

image_10m = image_30m.resample('bilinear').reproject({
'crs': crs,
'scale': 10.0})

image_10m_band2 = image_10m.select('B2')

get the projected image's projection and print

newproj = image_10m_band2.projection().getInfo()

print('New projection', newproj)

I also tried another version

transform_new = [
10.0,
proj['transform'][1],
proj['transform'][2],
proj['transform'][3],
-10.0,
proj['transform'][5]
]

projection_new = {
'crs': crs,
'transform': transform_new,
}

image_10m = image_30m.resample('bilinear').reproject(projection_new)

When I check the reprojected image's projection, I expected to obtain the projection in 10m, but only I got is a expection as follows:
ee.ee_exception.EEException Projection:
Argument 'crs': Invalid type. Expected: String. Actual: Type>.


texts: Cut layer into many parts using shapefile?

texts:
I have a layer which is a map of a town, and a shapefile which shows the boundaries within that town.

I wish to cut the layer up along each of the internal boundaries so that I end up with many layers, each containing the map of one small area as defined by shapefile. I know that I can use the clip tool to cut the map to the shape of the whole shape file, but I don't know how you use the boundaries contained within the shapefile to cut the map.

I am using ArcGIS Desktop 10 SP5 with the ArcView License


texts: Is it possible to create points that represent lines in QGIS 2.8.Wien?

texts: I want to create points to represent lines from a line layer. This is similar to creating points to represent polygons. Does anyone have an idea on any useful plugin or the best way to do this?


texts: How to place labels at the boundary of the current map extent

texts:
I have a layer of reference grid lines on top of several map layers. The grid lines should be labeled with coordinate values: vertical grid lines with x-values and horizontal lines with y-values. The label position should be at the border of the current mapcanvas extent.

I was unable to write an expression which uses mapCanvas().extent() due to illegal characters (the dot) . Geometry properties only will not return these coordinates.

My idea so far: creating point features at the intersection of grid lines and map extent triggered by extentChanged event, and labeling these points.

Has somebody a better idea?


texts: Connecting labels and points with arrows

texts: I found a previous thread on this and I am following this option.  I am hoping someone can split this into greater steps for me.  I am not certain of the detailed steps for creating a Generator Expression.  Can someone kindly expand on this process?


texts: Different data formats for the same graphic in Flex application

texts:
I have two ways to display an InfoWindow with the attributes of a graphic. The first is using a IdentifyTask to return the graphic and the other gets the selectedItem of a Datagrid (rovGrid) and using the findGraphicByAttribute code to return the graphic.

However, the format of the Date field is different in the returned graphic. The IdentifyTask graphic attribute shows "2/10/2005" and the Datagrid graphic attribute shows "1107993600000". Why are the two methods showing different date formats for the same feature?

Here's the partial code for each method

//IdentifyTask

identifyTask.execute(identifyParams, new AsyncResponder(onResult, onFault, clickGraphic));

function onResult(results:Array, clickGraphic:Graphic = null):void
{
    if (results && results.length > 0)
    {
        var result:IdentifyResult = results[0];
        var resultGraphic:Graphic = result.feature;

//Datagrid click

private function rovGrid_Click():void
{
    var obj:Object = rovGrid.selectedItem;

if (obj != null)
{
    MainMap.infoWindow.hide();

    var graphic:Graphic = findGraphicByAttribute(obj, layerROVGraphics);


texts: Spatial join and sum aggregation of two polygon GeoDataFrames in GeoPandas

texts:
I would like to join two GeoDataFrames with polygon geometries (purple and green, see screenshot below). However, I would like to get the sum of all intersecting green areas per purple feature they are intersecting with.

I tried
import geopandas as gpd

gpd.sjoin(purple, green, how="inner", op='intersects')

but that will just tell me which features are intersecting rather than returning the sum of the intersecting green areas. How would I do this in Python?


texts: Selecting through one-to-many relationship

texts:
I'm trying to make a species distribution project where I can select a species from a non-spatial list of species and the selection would highlight (select) all the counties to which the species has been attributed.

Simplified example data:
species, a non-spatial table of species (CSV with header):
"species"
"Mus borealis"
"Mus australis"

species_distributions, a non-spatial relation table (csv with header):
"species","county"
"Mus borealis","North"
"Mus borealis","Central"
"Mus australis","South"
"Mus australis","Central"

counties, a spatial feature class representing county boundaries (first row = field name):
"county"
"North"
"South"
"Central"
"West"

So if I select "Mus borealis" from the species table, the "North" and "Central" counties would be highlighted in the map. Seems like a simple scenario, but I have been unable to uncover a solution on my own. The closest I've found is a plugin "Select by relationship" but it doesn't appear to be available for QGIS 3 (I'm using QGIS 3.4.3). I have a PostGIS server available.


texts: Google Maps from QuickMapServices/OpenLayers Plugins won't show after zooming in and out

texts:
I use Google Maps (e.g. Google Satellite) as a basemap but after a few round of zooming in and out, the map would disappear. I have tried using maps from both OpenLayers and QuickMapServices (as the former seems to give problems to more people). QuickMapServices' maps were slightly better but the problem persists. I've checked that the CRS for the map layer and other layers match so don't think it's something to do with that.

Once the problem strikes, there seems to be nothing I can do. If I delete and add the map again, it will still not show. I can't even add other maps such as Bing Satellite. However, if I start a new project, the maps could work as per normal for the new project. This means that I have to abandon the whole project and start on a new one, and transfer whatever I've done to the new one. However, it happens so frequently that transferring is not viable.

Not sure if this will be useful but the version I'm using is QGIS 2.18 for mac.

I'm very new to QGIS.


texts: insert cursor "Row contains a bad value"

texts:
I'm having a problem copying data from a shapefile to a feature class using an insert cursor.  Here is the part that is causing trouble:

newRows = arcpy.InsertCursor(Target)
sourceRows = arcpy.SearchCursor(Christiansburg)

for source in sourceRows:
    neue = newRows.newRow()
    neue.setValue(TName, source.getValue(CName))
    newRows.insertRow(neue)

It takes data from a field in the source row, and then copies it to the new row.  The trouble is, when the source row has a forward slash in it, it fails when I try to do the insert, saying "The row contains a bad value.". If I edit the forward slash out of the source data, it works fine. If I replace the variable with a hardwired string that has a forward slash in it, it also works fine, so I am not 100% sure that it is the forward slash causing the problem.  I have tried using the str function on the result of getValue, and it still doesn't work. The field being copied from is street names, and the row causing trouble is a street that changes name.  If worst comes to worst, I could simply edit out the forward slash permanently, but I would rather not if I can help it.

Does anyone know what I am doing wrong here?

EDIT: Derped on copy/paste of code snippet. Should be fixed now.


texts: Summing a single column in attributes table QGIS 3.4

texts: I'm trying to sum up all of the values of a single column in QGIS3.4. I found some instruction on how to do this using the statistics panel, however, when I select the column I need to sum no values show up (see image). All the values in the column are numerical integers. I don't know if there's something wrong with my computer not playing nice with QGIS or if this tool is not working as it's supposed to (personally I'm assuming operator error).


texts: How to create a folder in HOME Folder with ModelBuilder in ArcGis 10?

texts:
I am just getting started with ModelBuilder and I would like to create a folder in the same folder where the respective mxd is located.
The tool have to be usable in different projects, so setting Current Workspace to a specific folder doesn't do the trick.

Is there any chance to get access to a path of the Home folder which would be the same in any given project?


texts: Difference between OGRGeometry and Geometry in C++ and python

texts:
According to the documentation, there is a python version of OGRGeometry's C++ class (http://www.gdal.org/ogr/classOGRGeometry.html) which is http://gdal.org/python/osgeo.ogr.Geometry-class.html

However, I don't understand why there is a lot of differences between these two class.

For instance, in Python there aren't any import method (importFromWkt).

Is it normal or maybe I'm on on the wrong way.

Thanks


texts: How do you remove the 1m by 1m gap polygons such as the one at the lower right corner?

texts:
How do I remove the 1m by 1m polygon in the lower right corner?

I am using ArcMap 10.4 and its Raster calculator.


texts:
I have coordinates of places that I would like to add to ArcMap. The coordinates are not written in long-lat format, so I am finding it difficult to use the "add xy" data command. For example, one coordinate is 7°35′13″N 1°56′06″W .

How can I add such coordinates to my shapefile?


texts: How to get Quantum GIS display a line or polygon feature created by Python package pyshp

texts:
I am not sure if I have discovered a regression, or if I have not enough sleep. While I manage to create shapefiles with point geometries with pyshp, line and polygon features don't show up in QGIS. In ArcMap, however, they show up as expected.

Steps to reproduce

Get pyshp (version 1.1.7)
Execute one or both of these example scripts, taken directly from the pyshp documentation:

Polygon

import shapefile
w = shapefile.Writer(shapefile.POLYGON)
w.poly(parts=[[[1,5],[5,5],[5,1],[3,3],[1,1]]])
w.field('FIRST_FLD','C','40')
w.field('SECOND_FLD','C','40')
w.record('First','Polygon')
w.save('test-polygon')

Polyline

import shapefile
w = shapefile.Writer(shapefile.POLYLINE)
w.line(parts=[[[1,5],[5,5],[5,1],[3,3],[1,1]]])
w.field('FIRST_FLD','C','40')
w.field('SECOND_FLD','C','40')
w.record(FIRST_FLD='First', SECOND_FLD='Line')
w.save('test-line')

Open the resulting shapefile in an empty QGIS session (tested in 1.7.3, 1.8.0 and 1.9.0 Alpha)

Expected result
A polygon or a line in the map window (at least after a "zoom to layer extent").

What I get
An empty map.

Did I miss something?

Update

After contacting the package author, version 1.1.9 was mentioned on the homepage. After an upgrade, the issue vanished. That's good customer support!


texts: Converting multiple drawing files to shapefiles

texts:
I have multiple dwg files and I need to convert them from CAD file to shapefile using Arcpy. I am trying to export each layer in .dwg file (i.e. annotation, point, polygon, polyline and Multipath layers, etc.) to shapefile as contains subsets of data within it.

I don't know how to access each layer of drawing and convert all drawing files to shapefiles. I have searched all existing but did not get sufficient information.

Any help?

Here is my code but I know this is not correct. Even I am not successful to access and print each layer in CAD. My ultimate goal is to convert multiple files dwg files to separate shapefiles using ArcPy.

import arcpy
from arcpy import env

env.workspace = r"C:DATA\DWG"
datasets = arcpy.ListDatasets()
for dataset in datasets:
    filesets = arcpy.ListFiles(dataset)
    for f in filesets:
            print f.name


texts:
Trying to get this code working with one OSM layer and one yahoo layer for OpenLayers 3 and Geoserver. Anyone have any idea where the problem migth be?

If I remove everything with the yahoo map, it all works fine (with only OSM map), however, i need the yahoo map to work as well.  Please help :)

P.s. src="http://api.maps.yahoo.com/ajaxymap?v=3.0&appid=euzuro-openlayers" script is loaded at the top, do not know if this might be the problem? Wrong script?

Open street Maps and Geoserver

<script>
//Defining projections
    var geographic = new OpenLayers.Projection("EPSG:4326");
    var mercator = new OpenLayers.Projection("EPSG:900913");

//Defining bounds
    var world = new OpenLayers.Bounds(-180, -89, 180, 89).transform(
        geographic, mercator
    );
//Defining map center
    var lund_center = new OpenLayers.LonLat(13.19, 55.71).transform(
        geographic, mercator
    );

    var options = {
        projection: mercator,
        displayProjection: geographic,
        units: "m",
        maxExtent: world,
        maxResolution: 156543.0399,

    };

//Defining main variables
    var map = new OpenLayers.Map("map", options, { controls: [] });

    var osm = new OpenLayers.Layer.OSM();
    var yLayer = new OpenLayers.Layer.Yahoo("Hybrid");

    map.addLayer(osm);
    map.addLayer(yLayer);

//List of layers

    var roads_wms = new OpenLayers.Layer.WMS(
        "Roads",
        "http://localhost:8080/geoserver/wms",
        {
            layers: "lu:vagar_alla",
            transparent: "true",
            format: "image/png"
        },
        {isBaseLayer: false, visibility: true}
    );
    map.addLayer(roads_wms);

    var raildroads_wms = new OpenLayers.Layer.WMS(
        "Railroads",
        "http://localhost:8080/geoserver/wms",
        {
            layers: "lu:railroads",
            transparent: "true",
            format: "image/png"
        },
        {isBaseLayer: false, visibility: true}
    );
    map.addLayer(raildroads_wms);

    var buildings_wms = new OpenLayers.Layer.WMS(
        "Official buildings",
        "http://localhost:8080/geoserver/wms",
        {
            layers: "lu:bygg_off",
            transparent: "true",
            format: "image/png"
        },
        {isBaseLayer: false, visibility: true}
    );
    map.addLayer(buildings_wms);

    var other_build_wms = new OpenLayers.Layer.WMS(
        "Other buildings",
        "http://localhost:8080/geoserver/wms",
        {
            layers: "lu:bygg_ovr",
            transparent: "true",
            format: "image/png"
        },
        {isBaseLayer: false, visibility: true}
    );
    map.addLayer(other_build_wms);

   ///////////////////////////////////////

//Map center and zoom
    map.setCenter(lund_center, 14);

//List of controls
    map.addControl(new OpenLayers.Control.LayerSwitcher());
    map.addControl(new OpenLayers.Control.KeyboardDefaults());
    map.addControl(new  OpenLayers.Control.Navigation());
    map.addControl(new OpenLayers.Control.ScaleLine());

//Coordinates of mouse position
    map.addControl(new OpenLayers.Control.MousePosition());

</script>


texts: QGIS - Use one layer to block another layers labels

texts:
I have two layers, one with points that I have set a label on and another with polygons which I would like to use to block the labeling of the points. According to this old answer to Avoiding labeling features if overlapped by another layer in QGIS? I could use the blocking of labeling option of the polygon layer, but this only works partially.
With "No labels":

With "Blocking":

Is this a bug or have I understood the option of "Blocking" wrong? (What other option than remove the points under the polygon do I have?)


texts:
I know that the shapefiles seen in figure 1 contain data taken from a database on the internet which is constantly changing.

However whenever I try change its properties such as its symbology, seen in figure 2 I am unable to change the symbols shape, colour, size etc

Why is this? Or is there a way to make this possible?

I use QGIS v3.2


texts: Remove the overlapping feature in ArcGIS

texts: I am using ArcGIS Desktop 10.7.1. I am looking for the removal of intersecting features. If any features intersect another feature then keep only one feature in the layer. In the following image, there are A-H features in a layer.  A & B and G & H are intersecting features. Among them, I need only one feature ( A or B and G or H) and all the remaining non-intersecting features.


texts: Calculating length and area with Google Maps API?

texts:
Google Maps API v2 had great methods for lines and polygons to calculate length and area.

v3 is lacking these methods (maybe just for now?).

What other options do I have for calculating line length and polygon area?


texts: How to deal with empty corners in raster when calculating fragmentation index

texts:
I have clipped buffer area around coordinate from one big raster (forest cover). Such created raster was loaded in R where PatchStat function (library SDMTools) gives fragmentation index.

After detailed look of such created raster I noticed that I have equal classification for non-forest area (0) and corners (also 0).

  My question is whether such surplus of empty area make an influence on the value
  of fragmentation index. I think solution bellow may be helpful, but I do not know how to do it.


texts: Finding lat/lng weighted midpoint from a list of points using QGIS

texts:
I am a truck driver and I have an historical list of cities, lat/lng, count where count is the number of times that I've made a delivery in that specific location.

I want to find the optimal location (lat,lng) where I should live, considering the history in this list. So that I'd be in a (statistically) good position for my usual work, should it not change too much in the future.

I've found this link, that explain the geographic weighted midpoint http://www.geomidpoint.com/calculation.html and it looks fine but maybe there are better approaches that I'm not aware of.

Is it a correct approach considering that I don't need road distances? (working on the sphere distance is good enough)
Is there an easy way to import such a CSV file into QGIS and find some kind of weighted  midpoint without relying on my terrible math skills?


texts: Categorizing most often attributes in column

texts:
I have a shp point layer with a "USER_ID" column. There is about 19 thousands features. Some of "ID"'s enters the list only ones or twice, but some of them enters 50 or more than 100 times.

I need to find the users which most often meet in this column. By all means it is possible to categorize this layer, but I think that the result shouldn't be clear. Is there a way to categorize only 50 most frequent values of a column?


texts: Perform Multiple Point in Polygon Queries PostGIS

texts:
I was wondering if there was a way to send one request to a PostgreSQL server with the PostGIS extension that contained multiple points, and then find the polygons that they are located in. For example, I am currently using the code below to query to see where multiple points are located in a polygon shapefile.

from sqlalchemy import create_engine, func
from sqlalchemy.orm import sessionmaker
from land_use_db_model import SpatialLandUse

# Query for unique lat lon vals
for i, df_tuple in enumerate(coord_df.itertuples(index=False)):

    lat = df_tuple[0]
    lon = df_tuple[1]

    """How to query to see where the point resides"""
    query = session.query(SpatialLandUse).filter(
        func.ST_Contains(SpatialLandUse.geom, 'SRID=4326;POINT({} {})'.format(lon, lat))
    )

    for result in query:
        coord_df.iloc[i, 2] = result.gridcode

So as you can see I am using sqlalchemy to perform the queries, but I am comfortable with using psycopg2 as well if that would be easier. I just feel like I spend a lot of time in this for loop and that it might be faster to send all of the points for the server to process at once.


texts:
I'm trying to convert a .gdb to .shp and I understand this is possible with 3.4. I didn't have any luck with 2.18.

I'm aware this question has come up previously but I have tried the other solutions and they haven't helped.

I am trying to get 3.4 standalone up and running and I continue to get the message/s

  The procedure entry point ... could not be located in the dynamic link
  library C:Program Files\QGIS\bin\qwt6.dll.

Further, three prompts appear asking for qscintilla2_qt5.dll, qca-qt5.dll and Qt5Sql.dll.

I have tried removing all previous traces of 2.18 using regedit. I have installed again as Administrator. I have removed virus software. I have tried copying the various .dll files over to the qgis.exe directory (btw, I can't see qgis.exe).


texts:
I imported a shapefile into QGIS, made a copy of it, and edited that copy. My edits consisted in changes in the coordinates of a few points. I then exported this edited shapefile as a .csv. For peace of mind, I imported the .csv to check my edits but, to my surprise, there were no changes although I saved the edits.
If I export the edited shapefile as a .shp everything is fine but as a .csv it doesn't seem to recognize the edits. As anyone come across this problem?


texts: pyqgis: issue with copying attributes of a features to a new created one

texts:
Context: the script selects features from a layer which overlaps another one above a certain percentage and creates a new layer with these selected features. I do not manage to copy the attributes from the attributes table, only the ID is reported.

What am i doing wrong?

Here is the code:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#for Severity
feats_lyr1 = [ feat for feat in layers[0].getFeatures() ]

#for Basins800
feats_lyr2 = [ feat for feat in layers[1].getFeatures() ]

selected_feats = []

print "processing..."

for i, feat1 in enumerate(feats_lyr1):
    for j, feat2 in enumerate(feats_lyr2):
        if feat1.geometry().intersects(feat2.geometry()):
            area1 = feat1.geometry().intersection(feat2.geometry()).area()
            area2 = feat2.geometry().area()
            crit =area1/area2
            if crit > 0.8:
                selected_feats.append(feat2)

epsg = layers[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "mem_layer",
                           "memory")

lyr = mem_layer.dataProvider()

for i, feat in enumerate(selected_feats):
    feat.setAttributes(selected_feats[i].attributes())
    lyr.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
print "done..."


texts:
I am working on google-maps kml. I have several shapefiles and using QGIS i have converted it to kml files. When i export the kml (save as kml in Qgis). It contains tag like <SimpleData name="Name">7</SimpleData> but i want tag structure like <Data name="name">
                        <value>7</value>
                    </Data> I have export the shape as kml and update the kml change style add BalloonStyle for that i need to modify structure otherwise it will not display info in popup window.Here is the link for rules of tag structure to display popup in google map as per kml extended data  If we use <SimpleData> tag then popup will not display lables in google map
When i change the structure of tags to <Data><Value>it will works gr8 ,Expected popup window is  I just want to change the tag structure.

Is there an automated way to change the tag structure from kml file?

For example, like below <simppledata> to <data><value>

<SchemaData schemaUrl="#new_circle">
                <SimpleData name="Name">7</SimpleData>
                <SimpleData name="gid">1</SimpleData>
                <SimpleData name="active_sub">4596</SimpleData>
                <SimpleData name="net_additi">5181</SimpleData>
                <SimpleData name="gross_addi">3917</SimpleData>
                <SimpleData name="churn">4617</SimpleData>
            </SchemaData>`

to

<Data name="name">
                    <value>7</value>
                </Data>
                <Data name="active_sub">
                    <value>4596</value>
                </Data>
                <Data name="net_additi">
                    <value>5181</value>
                </Data>
                <Data name="gross_additi">
                    <value>3917</value>
                </Data>
                <Data name="churn">
                    <value>4617</value>
                </Data>

here is the sample link in which i have commented the simpledata tag and add manually data and value tag.expected kml output To do it manually will take maximum time because i have thousands of records in other layers.


texts: I am working on a file which contains points from different areas of country. I have symbolized points. I want the feature count of respective area in a database form in .dbf file or excel etc. How can i get this count number (shown in figure) the dbf file.


texts:
I have QGIS an, access database and a Haicom gps all connected and running on a windows 8 machine.

I want to integrate all 3 items so that while in the field i can enter information into the database ie bird names, numbers etc AND automatically update XY location from the GPS ie a site AND see where i am in QGIS relative to other layers such as vegetation aerials etc.

status - at the moment QGIS easily logs a route and shows me where I am.

my question is - should i migrate access data to spatial database such as SpatiaLite or similar and then work on automatic logging of GPS data with a field entry

or write a macro/module in the access database to GET XY data from the GPS


texts: Calculate non overlapping areas with QGIS

texts:
I'm working with QGIS 3.16 Hannover under MacOS 10.13 environment. I've got two polygon layers: natural_parks (multiple features) and region (single feature). Some natural_parks features lay entrely within the region feature, some lay completely outside and some are only partly within. I want to calculate the areas of the natural_parks features that do not overlap with the region feature, assigning a 0 value for those enterely within. Actually, what I need is calculate the marine area (in hectares) for each natural park. I built the following code in the natural_parks layer with the field calculator, but it doesn't work
case when overlaps($geometry, geometry(get_feature('region', 'id', 1))) = 1
     then ($area - area(intersection($geometry, geometry(get_feature('region', 'id', 1)))))/10000
     when overlay_disjoint('region')
     then $area/10000
   else 0 end

What's wrong with my code?


texts: Using an asterisk in an annotation textstring

texts: I try to add an asterisk to the annotation (converted to shapefile), to indicate in my report that there is missing data for this polygon. However when I add the asterisk to the string and press "Toepassen" (in English it would be something like apply) nothing happens. How can I make the asterisk appear in the annotation?


texts: License Error when Reading or Writing to Esri Geodatabase in FME

texts:
After Upgrading  ArcGIS from 10.2.2 to 10.3, FME 2015 fails to run any workbench and shows this massage.

“Unable to perform licensing related tasks. Please ensure that ArcGIS is correctly installed
A fatal error has occurred. Check the logfile above for details
... Last line repeated 2 times ...
|ERROR |
A fatal error has occurred. Check the logfile above for details
|ERROR |
A fatal error has occurred. Check the logfile above for details”

The massage shows that the FME 2015 Do not recognize the ArcGIS 10.3 licence while my ArcGIS is working probably, I just found a link that talks about the same issue which is:

http://fmepedia.safe.com/articles/Error_Unexpected_Behavior/FME-cannot-read-or-write-to-Esri-Geodatabase-license-error-messages

and I just followed the instuction form ESRI about licensing ArcGIS product with Environment Variable but nothing helps

Have any one upgraded to ArcGIS 10.3 and faced the same problem ?


texts: Launching Processing tool User Interface using PyQGIS

texts:
I am trying to launch the UI of a processing algorithm using Python script but I cannot find a simple way for it (aside double clicking on it in the process toolbox of course). The plugin is a water network simulator and I just want to launch its pre-processing UI without using process.runalg() because this needs to enter arguments and I don't need it.

Some ways would be to trigger processes available in process toolbox or sending command to commander from python console. but I have no clue how and I am searching for 2 hours now, no results.


texts: Getting multiple functions to work with mouseover event in Leaflet

texts:
I am making a web mapping app using leaflet and some geoJSON data for a grid. I want to be able to mouseover each grid cell and both highlight that cell AND show a popup of the coordinates. In the code below, I have basically created two functions that should be triggered when the mouseover event occurs: highlightFeature(e) and coordinate(e). After my getJSON request to grab the geoJSON data for the grid, I am using layer.on({mouseover: <function call>}) to call the desired function when I mouseover a feature/grid cell. If I just use one function, it works fine. But if I try to do something likemouseover: highlightFeature, coordinate,` only the first function works. What am I doing wrong?

function highlightFeature(e) {
      var layer = e.target;

      layer.setStyle({
        weight: 0.5,
        fillColor: "#ffff00",
        dashArray: '',
        fillOpacity: 0.5
      });

      if (!L.Browser.ie && !L.Browser.opera) {
          layer.bringToFront();
      }
    };

function resetHighlight(e) {
  var layer = e.target;

  layer.setStyle({
    weight: 0.5,
    dashArray: '',
    fillOpacity: 0
  })
};

function zoomToFeature(e) {
  map.fitBounds(e.target.getBounds(15));
};

var popup = L.popup();
function coordinate(e) {
  var layer = e.target;

  popup
    .setLatLng(e.latlng)
    .setContent("Coordinates: " + feature.properties.LON + ", " +   feature.properties.LAT.toString())
  .openOn(map);
};
console.log(coordinate);

$.getJSON('./grid.geojson', function(data){
// add GeoJSON layer to the map once the file is loaded
  L.geoJson(data, {
    style: function(feature){
      return {
        color: "#000000",
        weight: 0,
        opacity: 0,
        fillOpacity: 0
      };
    },
    onEachFeature: function( feature, layer ){
      layer.on({
        mouseover: highlightFeature, coordinate,
        mouseout: resetHighlight,
        click: zoomToFeature
      });
    },
  }).addTo(map)
});


texts:
I want to run a LISA with a time series spatially lagged variable using PYSAL. I see that LISA will allow the definition of y & w but not Wy(t-k).

Any thoughts how I can get the lagged variable into PYSAL?

I also notice that PySal has a BV Morans but I can find nothing for a BV LISA.


texts: Python to Postgres insert Geometry syntax

texts:
Having trouble creating geom columns and inserting geometry lat lng values into postgres table

Is my syntax correct in creating the point geometry fields? st_makepoint is the point constructor. when I create the column should I be using that?
or should I be using the st_makepoint when inserting into the column?

*python 2.7

conn = psycopg2.connect("connected")
cur.execute('''CREATE TABLE example(pgeom geometry(POINT, 4326),
         dgeom geometry(POINT, 4326));''')
for x in csvfile:
    extracting info...
    1geom = lng1+lat1 #int values
    2geom = lng2+lng2 #int values
    query = '''insert into example(geom1, geom2) values (%s,%s);'''
    data = (1geom, 2geom)
    cur.execute(query, data)
    conn.commit()


texts:
I have a huge database with Sightings from whales and dolphins from 1995 till today.
So I want to sort the database with python.

Till now I can make FeatureLayers for each year

      for year in range(1995,2014):
      arcpy.MakeFeatureLayer_management
      ("Sightings95-14", "test%s" % year, "DATE >= date'01.01.%s' AND DATE < date'01.01.%s'" % (year, year+1))

now I want to sort it into each month (no matter of the yaer), because there are way more sightings in may than in december for example.

If I excecut this code:

      for month in range (01,12):
      arcpy.MakeFeatureLayer_management
      ("Sightings95-14", "test%s" % month, "DATE >= date'01.%s.1995' AND DATE < date'31.%s.2014'" % (month, month+1))

it makes the 12 FeatureLayers but in each Layer is the hole database included so he did not "sort" it right.

What is my mistake ?


texts: Characteristics of layer only for overlapping area of another layer

texts: I have got two layers, one with the land cover of Zimbabwe and another with the administrative divisions. I need to isolate a district and identify its land cover. I managed to cut out the district I need and make a new layer with it. However, every time I show the layer with the land cover, I get the land cover for the whole country. As, visually speaking, I only need that district, how can I do that? I tried to merge them or use "Union"/"Intersect", but all the land cover features disappear with these tools.


texts: How to convert Autocad lines to Esri polygon shapefiles using FME?

texts:
I am using FME to transform a dwg file into a polygon shapefile (.shp).

The dwg file is an urban plan of a city consisting of building blocks (type: polylines) and the properties within each block (type : line).

I used the GeometryCoercer transformer to convert type:line to type:polygon but it does not seem to work. I have also tried AreaBuilder, no success.

For those familiar with the FME software can you suggest Any solutions?

UPDATE

I attach an image of the dwg file (INPUT)

green closed polylines represent building blocks
orange lines represent properties

I used the workspace as provided by @Mapperz, but it didnt work out.
I set the geometry coercer  type : fme_line , then 2dForcer and AREA BUILDER writing to a .shp with attributed geometry shape_polygon. This is what I get.

When I set the output .shp geometry to shape_polyline I get the following

I need a shapefile where all buildings and all properties will have a polygon geometry in order to add attributes to each polygon.


texts: View or set seed value for Random point generator

texts:
I have a requirement to generate random points inside polygons but this generation must be repeatable.  Is there a way to either get the seed value that the random point generator is using or insert my own?

Got it.  works perfectly.  I was still running Essen.  Updated and works well.


texts: Is it possible to do drop markers on the map like in google maps (https://google-developers.appspot.com/maps/documentation/javascript/examples/marker-animations-iteration) in leaflet


texts: How to get longitude and latitude in degree units

texts:
When i navigate on my map , I receive the latitude and longitude, but I want them to be specified in degrees. I don't know what should I do please help
this is my code:

        function toggleControl(element) {
            for(key in drawControls) {
                var control = drawControls[key];
                if(element.value == key && element.checked) {
                    control.activate();
                } else {
                    control.deactivate();
                }
            }
        }

        function allowPan(element) {
            var stop = !element.checked;
            for(var key in drawControls) {
                drawControls[key].handler.stopDown = stop;
                drawControls[key].handler.stopUp = stop;
            }
        }

         function handleMeasurements(event) {
        var geometry = event.geometry;
        var units = event.units;
        var order = event.order;
        var measure = event.measure;
        var element = document.getElementById('output');
        var out = "";
        if(order == 1) {
            out += "measure: " + measure.toFixed(3) + " " + units;
            //alert("measure " + measure.toFixed(3) + " " +
             //                   + units + " ");
        } else {
            out += "measure: " + measure.toFixed(3) + " " + units + "<sup>2</" + "sup>";
        }
        element.innerHTML = out;
    }

    function toggleControl(element) {
        for(key in measureControls) {
            var control = measureControls[key];
            if(element.value == key && element.checked) {
                control.activate();
            } else {
                control.deactivate();
            }
        }
    }

    function toggleGeodesic(element) {
        for(key in measureControls) {
            var control = measureControls[key];
            control.geodesic = element.checked;
        }
    }

    function toggleImmediate(element) {
        for(key in measureControls) {
            var control = measureControls[key];
            control.setImmediate(element.checked);
        }
    }


texts: Spatial join points and polygons in PostGIS?

texts:
I have two spatial data sets in PostGIS tables:

polygon grid (table: grid_jakarta)
points (table:tweets)

For each point in 'tweets', I want to retrieve the reelvant attribute 'gridID' from the spatially corresponding polygon contained in 'grid_Jakarta'. Then I want to generate a new table that holds both. I assume this is a spatial join?

This is the query I have tried (based on a tutorial: https://datamodelprototype.wordpress.com/2014/01/27/linking-point-data-to-polygon-data-using-qgis-or-postgis-sql/), but to no avail:

SELECT pts.*, grid.gid as gridID
FROM "tweets" AS pts
INNER JOIN "grid_jakarta" AS grid
ON st_within(pts.location, grid.the_geom) AS result

Note, I need to use PostGIS due to the large quantity of points and polygons.

I have also checked previous QAs (e.g How to perform a spatial join of point and polygon layers in PostGIS?), but these aren't trying to achieve exactly what I'm trying to achieve.


texts:
I am using OpenLayers to render a WMS inside of a div tag.  This is right out of the OpenLayers documentation.  My WMS layer is rendering well.

Using the jQuery UI documentation I have been able to create a draggable div tag that will, hopefully, house controls for editing vector features.  Using the code below I am able to draw my draggable box.

JavaScript:

$('#makeMeDraggable').draggable();

HTML:

I was able to get overlap working by setting the map position to absolute and assigning appropriate z-indices to each div. (In the CSS file)

#map{ position:absolute; z-index:1;}
#makeMeDraggable{z-index:2;}

This allows the box to be dragged around the map's "container" (div tag) without issue.  I am struggling to figure out how to have the map's position remain constant as I drag the box?  The map is registering an onClick event over the entire page and so the box and map move in tandem.  I am relatively new to jQuery and OpenLayers.  I am probably missing something straight forward in the documentation?

Edit: URL to jsfiddle.net for anyone interested in checking out this problem.  The map displays after zooming once.  I still need to modify that...

http://jsfiddle.net/Dn3fw/


texts: Sum/aggregation of polygon's values within other polygons using dyplr


texts: Changing the spatial resolution - stacked raster in R

texts:
Want to change higher spatial resolution UAV dataset stacked in R. Used raster package to stack the individual band. Want to change the spatial resolution to 1 m.

stacked_UAV  <- stack(stack_list)
### Check the original resolution
res(stacked_UAV)
0.0289 0.0289
## Number of rows original stacked image
ncol(rast) <- stacked_UAV@ncols
ncol(rast)
5193
## Number of rows original stacked image
nrow(rast) <- stacked_UAV@nrows
nrow(rast)
13680

Tried following options using "raster package' but failed to get the target 1 m spatial resolution.

test_with_aggregate <- aggregate(stacked_UAV, fact=4)
res(test_with_aggregate)
0.1156 0.1156

Then tried to use the 'resample' function in the raster package but could not get the target spatial resolution due an error in setting spatial resolution to 1 m, this option is also not working.

### Empty stack
rast <- stack()
### Same spatial extend as the original stack
extent(rast) <- stack(stacked_UAV)
## Set resolution
res(rast) <- 1
Error in x@layers[[i]] : subscript out of bounds
stacked_UAV_1 <- raster::resample(stacked_UAV, rast, method='bilinear')


texts:
I'm looking for a roughly equivalent tool in QGIS to the Repair Geometry Tool in ArcGIS Desktop. When I run some common vector processing tools in QGIS, I occasionally run into invalid geometry errors and have yet to find a quick/easy fix that doesn't involve manual editing.

Is there an automated process that I can use?


texts: QGIS Time manager Changing Color Polygon

texts:
I'm using QGIS 3.6.3. Noosa.

I have 10 polygons, each has a starting time and an end time. Getting Time Manager to run was no problem at all. However, there are two things I would like to add in my dynamic visualization.

With each ongoing time frame I would like the color of the polygon to change, meaning that when in January 2019 the polygon is dark red, it should be less dark in February 2019 and by the end date in December 2019 be a very light red color.
Polygons disappear at the end of the time span, I would like them to be visible (as stated in 1, in a different color than at the beginning of the time span).

Does anyone have any suggestions?


texts: Rule-based color in dataPlotly bar chart in Atlas

texts:
I have created an Atlas in QGIS 3.10 and have made a bar chart that changes accordingly to the corresponding Atlas feature.

x-axis: area type and y-axis: area in square meters

I also managed to set the color for each bar to correspond to the color of the area type shown in the current atlas map by using an data defined expression for the bar color

case
  when "artype"=11 then color_rgb(255,85,255)
  when "artype"=12 then color_rgb(197,197,197)
  ...
end

Now, I wonder if it is possible to

change the label sizes (or colors and font types) without using
HTML-tag like  and  e.g.
Arealtype. I've tried
 and

without any luck.
just highlight the largest area type (the highest bar), i.e. showing
the largest area type in orange color, the rest of the bars in grey.
I have tried different variants of

case
   when attributes(get_feature('Dissolved','areal_u32',maximum("areal_u32",filter:="rn"=@atlas_featurenumber)))['artype']=11
then color_rgb(255,127,0)
   ...
   else color_rgb(229,229,229)
 end

but always (of course?) ending up with all bars being orange.
keep all area types on the x-axis in a predefined order. The
bar height will be zero for area types not present in the current
atlas map. This would make it even easier to compare the bar charts
for each atlas map.

Can all/some of this already be achieved or do I have to wait for QGIS 3.11?


texts: How to Convert db.ListProperty(GeoPt) to GeoJSON Polygon

texts:
My Python Google App Engine (no tag?) app has a datastore model of an Area of Interest, simplified as

import geojson
class AreaOfInterest(db.Model):
    name = db.StringProperty(required=True)
    #Geometry of area boundary
    coordinates = db.ListProperty(db.GeoPt, default=None) #convex polygon without holes

I wish to export this as GeoJSON, but am unsure of a (preferably Pythonic) way to convert the List of GeoPt to a geojson.Polygon object. Although I can convert to an array of Points with

coords = []
for c in self.coordinates:
    p = geojson.Point((c.lat, c.lon))
    coords.append(p)

I wish to embedd a Polygon into a GeoJson FeatureCollection thus:

def geojsonArea(self):
    polygon = geojson.Polygon((self.coordinates))
    #ValueError: Expect __geo_interface__' or a mapping, got datastore_types.GeoPt
    geojson_obj =   {
        "type": "FeatureCollection",
        "area_name" :self.name,
        "features": [
          { "type": "Feature",
             "geometry": {
               "type": "Polygon",
               "coordinates": polygon
             },
             "properties": {
                "name": "boundary"
           }
        }
        ]
    }
    return geojson.dumps(geojson_obj)

But the call to geojson.Polygon expects something like
 Polygon([[(2.38, 57.322), (23.194, -20.28), (-120.43, 19.15), (2.38, 57.322)]])

What is the cleanest way to make this conversion?
Is it possible without looping through the list while building a string to pass to Polygon()?


texts: Spatial relation (correlation) between polylines

texts:
I'm evaluating the spatial relationship between two polylines (morpholineament and magnetic lineament). If this two polylines are spatially correlated or not.
Traditionally, this analysis is realized by visual inspection:

The polyline A is close to polyline B, so they are spatially correlated

This is a subset of my polylines. In a visual inspection, the E-W and NE-SW sets are spatially correlated to the magnetic lineaments. However, this structures may shift, be dislocated or rotated, so this traditional method may not work properly.
Carranza (2009), for instance, uses the cumulative frequence of the distance between the mineral deposits (point) and fault system (polyline) to analyze the spatial correlation.

Carranza, E.J.M., 2009, Controls on mineral deposit occurrence
inferred from analysis of their spatial pattern and spatial
association with geological features: Ore Geology Reviews, v. 35, p.
383–400, doi:10.1016/j.oregeorev.2009.01.001.

I can "convert" one of my polyline to points and apply the Carranza's method, however I don't think this will be correct.
So, I would like to avoid the visual inspection and apply some statistic method, do you know any method to do so?


texts:
I've searched both the 2008 spec and RFC7946 and can't find anything about when two GeoJSON objects would be considered equal.  I wonder if anyone can enlighten me?

It seems to me the options are:

Two GeoJSON objects that contain the same geometry are never equal,
and are both distinctly unique features.
Two GeoJSON objects are equivalent if their geometry coordinates are
exactly equal, irrespective of their geometry type.
Two GeoJSON objects are equivalent only if their geometry -
both coordinates and type are exactly equal.

Option 1 poses questions on handling assertions about the same geometry from multiple sources, as well as uniquely identifying them.

Option 2 poses questions regarding handling multiple geometry type declarations, which effectively becomes a multi-valued property

Option 3 poses questions about a GeometryCollection asserting multiple geometries that contain the same coordinates but different types

The GeoJSON-LD proposals appears to suggest a combination of all three options are possible via @id, which means equality is outside of the RFC spec and not defined without an external property.

I have been encoding all my geometries as unique base64 data:urls with the polyline algorithm which works well, and have allowed geometry type to be multi-valued and merged all other properties.

My thinking is both the properties object and any foreign members are explicitly outside of the spec so can be safely merged, though I'm not sure how correct the semantics of that are.

Example below:

[{
    "type": "Feature",
    "geometry": {
        "type": "LineString",
        "coordinates": [[123,456],[456,789],[789,012],[012,234]]
    }
    "properties": {
        "name": "Bob",
    },
},
    "type": "Feature",
    "geometry": {
        "type": "LineString",
        "coordinates": [[123,456],[456,789],[789,012],[012,234]]
    }
    "properties": {
        "name": "Jane",
    },
},
    "type": "Feature",
    "geometry": {
        "type": "MultiPoint",
        "coordinates": [[123,456],[456,789],[789,012],[012,234]]
    }
    "properties": {
        "name": "Fred",
    },
},
    "type": "GeometryCollection",
    "geometries" : {
        {
            "type": "Polygon",
            "coordinates": [[123,456],[456,789],[789,012],[012,234]]
        },
        {
            "type": "MultiPoint",
            "coordinates": [[123,456],[456,789],[789,012],[012,234]]
        }
    }
    "properties": {
        "name": "Mary",
    },
}]

I don't have the rep to add tags for JSON-LD or GeoJSON-LD, perhaps someone else could for me.


texts: ArcGIS Python TableSelect_analysis ERROR 000210: Cannot create output

texts:
I have a very strange issue:
I try to execute the following script which uses TableSelect_analysis and construct as outputs the corresponding tables from another table depending on a specific attribute. For example I have 1000 rows with urban, rural and industrial areas and I want to output 3 files each with the corresponding kind of area.

The code:

for idx, description in enumerate(greeklishDescrList):
    description = description.replace(" ", "_")
    outSHP = outDir + description + u".shp"
    idxOriginal = descrList[idx]
    arcpy.TableSelect_analysis(r"C:\r*******d\R********l.gdb\l*********I.shp", outSHP, "\"descr\" = '" + idxOriginal + "'")

What happens is this:

In the first iteration of the for loop it does select all the rows with the correct value "descr".

1st ISSUE:

Output name in the gdb is "shp" when it should have named it as I have coded elsewhere as RURAL.shp . In the contents of "shp" table everything is as expected.

2nd ISSUE:

In the next iteration it breaks down showing the following message:

  ExecuteError: ERROR 000210: Cannot create output
  C:\r*********d\R********l.gdb\URBAN.shp Failed to execute
  (TableSelect).

Through the Python scripting environment of ArcGIS (I code through PyCharm) it also shows the following

Troubleshooting:

The path is correct. It is clear also from the fact that it constructes the problematic "shp.shp" table with -however- the correct values.
There are permissions available for the same reason as (1)
I have excluded the gdb and the whole directory from my Microsoft Security Essentials
The files do not exist (I believe!). I delete them from the gdb through ArcGIS and I restarted everything multiple times

ArcGIS 10.3.1 with Python 2.7.8


texts: GDAL CSV Driver with "GEOMETRY=AS_XYZ" option not creating X Y columns

texts:
I have an iOS GIS app written in Swift that relies heavily on GDAL 2.2 for the spatial data management.  I've written Swift wrappers around the parts of the standard-C GDAL that I want to use.   I'm adding options to the Swift wrappers to export data in various other formats.
However, when exporting to CSV, I'd like to include the coordinates (for point features).  However, the X/Y columns never get created, despite passing the "GEOMETRY=AS_XYZ" option to GDALCreate.  (I was also unable to find any examples when doing a web search - just lots of examples of doing it with the ogr2ogr utility.)
I've not used options with GDALCreate before, so it's possible that I'm doing it wrong.  There are also some complications with doing this from Swift.  In particular, converting from Swift string objects to standard C strings (char pointer) and especially from a Swift array of strings to standard C array of strings (pointer to char pointer).
Here's how I'm attempting to do this:
let gdalDriver = GDALGetDriverByName("CSV".cString(using: .utf8))
let options = ["GEOMETRY=AS_XYZ", "CREATE_CSVT=YES"]
var cOpts = options.map { strdup($0) }
cOpts.append(nil)
ogrDS = GDALCreate(gdalDriver, url.path.cString(using: .utf8), 0, 0, 0, GDT_Unknown, &cOpts)
for ptr in cOpts { free(ptr) }

From this point, I can copy GDAL features to it from my source points dataset, and that works fine (all features with all attributes appear in the CSV) except that there are no coorinates (and no columns for the coordinates).
Has anybody here had any success creating CSV data sets with GDAL standard C (or C++) including X/Y columns for points?  What am I doing wrong here?
(Should I post this in Stack Overflow instead of or as well as here?)


texts: Trying to convert tiff+tfw to mbtiles. Projection issues?

texts:
I have quite a number of huge TIFF files, each TIFF file has a TFW file associated. gdalinfo shows that the TIFF files are not real GeoTIFF (no embedded data). I know the TIFF's are in WGS84, UTM format. I was able some time ago to load these files into CompeGPS, which created properly calibrated images in that software.

Now I need them converted to MBTILES to use with AlpineQuest on Android.
I can create the mbtiles with this command:

gdal_translate test.tif test.mbtiles -of MBTILES

The problem is that the map is not correctly calibrated!

It should sit around 45N/7E, but if I open the MBTILES in QGis or in AlpineQuest I see it around 40N/1E.

I suspect the problem is with the source image, since also the TIFF+TFW will create a similarly misplaced image in QGis.

Now I am quite an amateur in GIS matters, so I am puzzled and really don't understand what I am doing wrong or what I should be actually doing?

Can somebody help me undestand what I am doing wrong?


texts: Need to create custom map based the spatial Database


texts: How to correctly display SVG path graphic after map scale change using ESRI JavaScript API?

texts:
I've built a web map that shows the locations of fleet vehicles based on their GPS data. I managed to add a nice feature that animates the routes these vehicles take on their way to destinations using SVG paths generated from street network data, however there is an issue with how the route graphics display after panning or zooming the map. Route calculations happen on my server-side app and are broadcast to my client app each time the vehicle updates its location. The first time a route line is drawn it looks correct for the current map scale, but if I zoom in the route line will look "clipped" or broken.

Here is an example of a route after the initial draw and animation. Imagine the yellow line originating from the black "5" graphic and snaking over to the red cross icon.

Here is the same route graphic a moment later zoomed in one level:

If I were to stay at that scale the line would never fully draw. Only after quickly zooming out and then in does the graphic "fix" itself and fully display. This issue is present whether I have SVG path animation on or not.

I believe this has something to do with how the SVG viewport is initially set and then not updated on a map scale change. One obvious quirk is that if a route graphic should be shown but it happens to be off the visible portion of the map it will not be visible even when the map is panned to that route location. Only after zooming in and out will the graphic finally render. I hope someone can offer a solution for making sure the route lines are refreshed and visible after scale changes. I will post some of the graphics code and what I've already tried. For the record I'm using ESRI JavaScript api v3.14 with Dojo v1.10.4 and testing on Chrome. My path animations were originally set only in CSS but I have tried Snap.svg for comparison and the added benefit of IE11 SVG support.

//Initialize ESRI graphics layer
routeLayer = new GraphicsLayer();

//My "RouteMaker" script handles getting new coordinates and creating the graphic...
var coordinates = routeObj.coords;
var routeLine = new Polyline(coordinates);
var sls = new SimpleLineSymbol('STYLE_SOLID', routeColor, 5);
this.graphicRoute = new Graphic(routeLine, sls);

//...which is then added to the graphics layer
var newRoute = new RouteMaker(routeObj);
routeLayer.add(newRoute.graphicRoute);

//Prior to rendering the route graphic I use the 'graphic-node-add' event to add
//shadow filter and svg animation. Here I am using Snap.svg to animate the line
var filter = filters.createFilter({x:"0", y:"0", width:"120%", height:"120%"},
    [
        filters.feGaussianBlur({"in":"SourceAlpha","stdDeviation":4}),
        filters.feOffset({"dx":5,"dy":5,"result":"offsetBlur"}),
        filters.feMerge("offsetBlur","SourceGraphic")
    ]);
routeLayer.on('graphic-node-add', function (evt) {
    var shape = evt.graphic.getShape();
    shape.setFilter(filter);

    var pathNode = evt.graphic.getNode();
    var snapPath = Snap(pathNode);
    var lenPath = snapPath.getTotalLength();
    snapPath.attr({
        'stroke-dasharray': lenPath + " " + lenPath,
        'stroke-dashoffset': lenPath
    }).animate({
        'stroke-dashoffset': 0
    }, 500);
});
//Alternatively I have tried leaving out Snap.svg and just using CSS, but the
//issue remains.
path {
    stroke-dasharray: 2500;
    stroke-dashoffset: 2500;
    animation: dash 5s linear forwards;
}
@keyframes dash {
    to {
        stroke-dashoffset: 0;
    }
}
//My attempt to force the route graphic to redraw. This seems to have no effect.
this.map.on('extent-change', function (evt) {
    console.log('map extent evt: ', evt);
    routeLayer.redraw();
});


texts: mosaicing geotiffs with gdal creates gaps due to alpha layers

texts:
I have a series of geotiffs which I'm trying to mosaic together, however when they are combined their individual alpha layers overlap sections of the other images and create gaps in the final mosaic:

My approach is first to generate a .vrt from all the geotiffs:

gdalbuildvrt -hidenodata mosaic.vrt geo_pict20140910_131*

And then to convert that .vrt to a geotiff with gdal_translate:

gdal_translate mosaic.vrt mosaic.tif

Is there a way I can mosaic only the non-alpha portions of these geotiffs to avoid these gaps?

Here's a link to the original geotiffs:  https://www.dropbox.com/s/udjq9bter9fmgxz/gtiffs.zip?dl=0

Here's a link to the .vrt and final .tif mosaics:  https://www.dropbox.com/s/oov9vatz69f31zg/mosaic.zip?dl=0

Here's the gdalinfo for the final mosaic:

Driver: GTiff/GeoTIFF
Files: mosaic.tif
Size is 4782, 2837
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-121.460884208999985,38.583432355500015)
Pixel Size = (0.000000591354474,-0.000000591354474)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-121.4608842,  38.5834324) (121d27'39.18"W, 38d35' 0.36"N)
Lower Left  (-121.4608842,  38.5817547) (121d27'39.18"W, 38d34'54.32"N)
Upper Right (-121.4580564,  38.5834324) (121d27'29.00"W, 38d35' 0.36"N)
Lower Right (-121.4580564,  38.5817547) (121d27'29.00"W, 38d34'54.32"N)
Center      (-121.4594703,  38.5825935) (121d27'34.09"W, 38d34'57.34"N)
Band 1 Block=4782x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=4782x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=4782x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=4782x1 Type=Byte, ColorInterp=Alpha


texts: Fetch the selected value from scale toolbar

texts: can somebody guide me how to fetch that value shown in the image? using jquery and or PHP language?


texts:
I can't understand the meaning of feature in QGIS QgsFeature class. How it is related to layer and geometry?

I have read QGIS documentation but I can't understand what is feature.


texts: Using REPLACE function in ArcGIS Server Query - Is it supported?

texts:
I'm working on a Flex application for a client that will be the replacement for a legacy application they have that was written in .net/adf. In their current code base, they are using an ADF query with a Replace function inside the where clause.

adfQueryFilter.WhereClause = "REPLACE(TMS,'-','') like '%" + ColValue + "%'";

The reason for this is they have multiple parcel tables, and each table has a TMS column that they are querying on. The TMS values for each table have unique number groupings such as ###-##-##-### or ####-##-##-## etc. The consistency is that they all contain 10 digits. The ColValue variable above will contain a sequence of numbers WITHOUT hyphens in it (example: ##########). ColValue only exists in their data this way and cannot be converted to a hyphenated format without configuration rules and coding. That extra config and coding is what I'm trying to avoid if possible, and the Replace function seems to be capable of doing just that - in ADF anyway.

So, the question then is, can REPLACE be used in an ArcGIS Server Query, or is it only supported by ADF?

The following query will generate a selection in ArcMap: Select * from ParcelTable where: REPLACE(TMS,'-','') LIKE '%63%'

However, when I try that same query on the ArcGIS Server rest endpoint it fails.

This is a ArcGIS 10.2.1 map service fyi.

Thanks for any input on this!


texts: Finding average distance from polygon to multiple points

texts:
I am trying to find the average distance from a polygon to multiple points on a map. Specifically from a parcel to "downtown hubs". I have used the "Near" tool to generate a table showing the distance from each parcel to each point. Now I want to average the distances for each point. Then join the result for each FID back to another feature class in a new field.

My current table is below

I want this to become something like:

INPUT_FID | DISTANCE

1  |   xxxxx

2 | xxxxx

3 | xxxxx


texts: Making PostGIS extension persist in the search_path of a PostgreSQL database?

texts:
Context

I have a script to set up a PostgreSQL database along with the postgis extension which I decided to install in an other schema than public, namely, the postgis schema. In my script, this schema is then set in the search path as explained here;

SET search_path TO "$user", postgis, topology, public;

After having successfully executed my creation script, I was totally surprised by this error while connected with the user role I defined (different than the postgres role) and using a PostGIS function to test it was correctly installed:

$ psql -U user -d my_database
psql (10.12 (Ubuntu 10.12-0ubuntu0.18.04.1))
Type "help" for help.

my_database=# SELECT st_astext(ST_SetSRID(ST_MakePoint(12.0, 42.0),4326));
ERROR:  function st_makepoint(numeric, numeric) does not exist
LINE 1: SELECT st_astext(ST_SetSRID(ST_MakePoint(12.0, 42.0),4326));
                                    ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Then I checked for the search_path;

my_database=# show search_path;
        search_path
---------------------------
 "$user", public, topology
(1 row)

Astonishing; PostGIS is no more here!

I have to enter the SET command again, manually, to make it come back:

$ psql -U user -d my_database
my_database=# SET search_path TO "$user", postgis, topology, public;
SET
my_database=# show search_path;
                 search_path
---------------------------------------------
 "$user", postgis, topology, public
(1 row)

But this doesn't persist across re-connections.

Question

How to make it persistent?

Related questions and links

Use PostGIS functions from an other schema (than public)

PostGIS functions do not exist even when public is in the search path

https://intellij-support.jetbrains.com/hc/en-us/community/posts/208393425-type-geometry-does-not-exist


texts:
I want to create Military Symbols in Flex Viewer 2. How I can accomplish it, including all the names and their code in an XML file and calling it within flex viewer is a pain. Is there any SDK available for this purpose. That can generate symbols for me with like MOLE provides.
If the SDK is available for JAVA or .NET, Please point it out, I will get an idea from it. I had heard about Defense Solution SDK and ARC Objects but want a bit guidance.


texts: Formatting longitude/latitude in arabic language?

texts:
How should longitude and latitude be properly formatted in the arabic language?

e.g. how should the following coordinates be displayed in arabic language:

40°26′46.23″ N, 79°56′56.12″ W?


texts: Right to left languages in GeoNetwork v3

texts: Is there any support for Right-To-Left (RTL) languages such as Farsi and Arabic in GeoNetwork v3?If yes, any document for implementing it?


texts: How to convert contours and raster surface from Petrel to shapefiles and rasters that can display in ArcGIS?

texts:
We have Petrel which is used to produce a surface of a geological horizon below the earth's surface. In Petrel we have contours (vectors) and a raster image of the same area, also visualising elevation of the surface.

Has anyone got any experience in ways to convert this vector dataset and raster dataset in Petrel to a format suitable for ArcGIS?

A previous employee sent it to Petrosys from Petrel, which could then output shapefiles of the contours and a .BIL raster, which I can display in Arc.

This, however, is something no current employee has experience with. FME does not handle the regular Petrel output file types of


texts:
I have this big land use map. I extracted shapes associated with mountains so that I can improve the classification (for example, I corrected clouds in the mountains to be forest, among other things). Now I want to flatten the two layers with the extracted layer over-writing the info on the original map. I tried union, merge and other things, but the info I want overwritten is still there (ex. the clouds are still there).

I checked for instance the result of the union and I saw that the all the data columns of both shapefiles were appended together. The effect I want is for the data of the top layer to replace the data of the lower layer. In Photoshop, this is just basic flattening (i.e. pixel value of higher layer replaces all pixel values of underlying layers).


texts: ogr2ogr: issue with batch clipping shapefiles with converted kml file as clip source

texts:
I wrote this script that (running in MSYS) converts a kml to shp and applies a re-projection to it. The resulting shp then has the same SRS as the shps that I want to clip with it. I use a for loop to look up the shp files in a certain folder and then want to apply the clipping.
The conversion throws some neglectable (?) warnings and yields a working shp-file.

Warning 6: Normalized/laundered field name: 'description' to 'descriptio'
Warning 6: Field timestamp create as date field, though DateTime requested.
Warning 6: Field begin create as date field, though DateTime requested.
Warning 6: Field end create as date field, though DateTime requested.
Warning 6: Normalized/laundered field name: 'altitudeMode' to 'altitudeMo'

However, the second part is not working (no errors, just the ogr2ogr usage infos as it appears if there are syntactical errors)

cd D:/GIS_DataBase/GIS_Tirol/Scharnitz

outdir=D:/GIS_DataBase/GIS_Tirol/Scharnitz/Data
indir=D:/GIS_DataBase/GIS_Tirol/Tirol_Verbreitungskarten/Verbreitungs_Daten

mkdir outdir

# reproject and convert clipsource to desired SRS and shpfile-format
ogr2ogr -overwrite -skipfailures -append -f "ESRI Shapefile" -t_srs EPSG:31254 clpsrcfile.shp scharnitz_gem_grenze.kml

# clip all files in indir and save to outdir
for x in $indir/*.shp
  do ogrinfo -al $x | grep Extent
  ls $x
  ogr2ogr -clpsrc clpsrcfile.shp $outdir $x
done


texts:
I have a problem with a shapefile (EPSG:31287 MGI_Austria_Lambert). I tried to solve it on my own, but I just have no clue as to what's the error.

So I use a basemap from basemap.at (WMTS: https://www.basemap.at/wmts/1.0.0/WMTSCapabilities.xml -> EPSG:3857 - WGS 84 / Pseudo-Mercator). When importing the basemap and then the shapefile in QGIS it looks like a weird distortion (Austria shapefile), the Polygons (which represent houses) are not exactly over the houses in the basemap, but a few meters away. When I do the same in ArcMap, there is no Problem and every Polygon matches the house in the basemap.

But I need to export it as KMZ or KML and when I do so, and then import it in Google Earth I get the same distortion like in QGIS)

I really hope someone has a clue or knows where the problem stems from.

ArcMap:

QGIS:


texts:
The Esri docs for PostgreSQL data types supported in ArcGIS v10.6 says in part [emphasis added]:

  If you create your geodatabase using the Create Enterprise
  Geodatabase geoprocessing tool, your database is created using a
  template other than the PostGIS template. Therefore, if you want to
  use the PostGIS geometry or geography type, you must create a database
  manually using the PostGIS template. Then you can specify your
  existing database when you run the Create Enterprise
  Geodatabase geoprocessing tool, and the geodatabase will be created in
  this PostGIS-enabled database.

How exactly does one "create a db using the postgis template"?

We intend to use the postgres GEOGRAPHY type for GEOMETRY_STORAGE instead of the Esri default ST_GEOMETRY (ref).

I've installed Postgres 9.6 using the Esri supplied installer* and used Stack Builder to install PostGIS extension. I'm able to connect to Postgres with pgAdmin from a remote computer.

* PostgreSQL_DBMS_for_Windows_963_158796

update: fixed installed postgres version typos. Added ArcGIS version.


texts:
I am looking for advice on how to set up a multivariate regression using raster data. My dependent variable is the location of firms in a specific industry.
My independent variables are:

lease rates
Within 1 mile of large firm
within 1 mile of incubator
within 1 mile of university
within 0.5 miles of freeway
In industrial zoned area

Using the band collection statistics tool, I can get correlations, but this does not give me the multivariate regression model I am looking for. I could export and run in another statistical software, but I can't figure out how to make the grid identical, and also the attribute table for rasters does not give coordinates or specify each cell so that I can match records across variables.

Any advice on this? I am open to any solution!


texts: Is there a maximum number of feature classes for an Esri File GeoDatabase?

texts:
Is there a maximum number of feature classes for an Esri File GeoDatabase ? I'm running a script that creates feature classes, however, after around 1000 I receive a schema lock error:

ERROR 000210: Cannot create output.


texts: How to join a WFS layer to a stand-alone table in OpenLayers

texts:
I'm a complete novice at OpenLayers, and am checking whether the following is possible.

I envisage a dataset of polygons in GeoServer, being served as a WFS layer into OpenLayers. The polygons won't have any attributes other than a unique identifier (eg State name).

The attributes for the polygons will be provided by a stand-alone table, which should match to the polygons via the unique identifier.

The join needs to be configured within the user's browser (rather than on the server) since the table will be unique to that user and that session - the table is generated dynamically as the result of the user's actions.

So my questions are:

is it possible to join a WFS layer to a stand-alone table?
is there any sample code showing how to set this up?

Thanks


texts: Linking excel/csv file to QGIS to update fields automatically using Macros

texts:
I am creating animated maps related to a couple of public health projects on Dengue and Zika. All our data is stored in Commcare, and then we can download it as an excel file, with all the data and headers (column titles).
So, if we need to create maps at any point, we have to simply export the data from commcare into excel, then convert that into a csv and join it to a shapefile of Iquitos called 'Noviembre_P01' (which contains the 'LOC_CODE' field on which is based the join with the 'loc_code' field in the csv file). The csv file is called 'project_2_gis_mapping'.

I want to automate this process of the attribute tables in qgis getting updated everytime I export this data from Commcare into excel as a csv file.

In the entire csv file, I have one column/field of interest - 'first_bloo'. I have created 1 additional column in my attribute table with this same name.

From my understanding, it seems that writing a Macros code will help me achieve this. Please can someone help with the Macros code for the same? Or, is there an easier/more efficient way of achieving this process?


texts:
having trouble with debuging QyPython plugins from Eclipse:
- Win 7 x64
- QGIS 2.7
- Eclipse Luna with PyDev installed

Tried to debug using Remote Debugger (http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/ide_debugging.html)
Remote Debbuer displays "Debug Connection Failed"
Tried to debug using settrace() as it is explained here:
http://linfiniti.com/2011/12/remote-debugging-qgis-python-plugins-with-pydev/
Few days this approach working fine but today when I started QGIS desktop -
desktop hangs on Startup Popup with "Restoring loaded plugin" status.
I.e. the break point in my plugn is not hit.
If I remove settracte() from code - QGIS is started OK.
Looks like some break in link between QGIS and PyDev Debug Server.

I'm using Simple Test Plugin created by plugin builder.
Here is the code for setrace()

import sys
sys.path.append("C:/eclipse/plugins/org.python.pydev_3.9.2.201502050007/pysrc")
from pydevd import *

FORM_CLASS, _ = uic.loadUiType(os.path.join(
    os.path.dirname(__file__), 'Test_dialog_base.ui'))

class TestDialog(QtGui.QDialog, FORM_CLASS):
    def __init__(self, parent=None):
        """Constructor."""
        super(TestDialog, self).__init__(parent)
        # Set up the user interface from Designer.
        # After setupUI you can access any designer object by doing
        # self.<objectname>, and you can use autoconnect slots - see
        # http://qt-project.org/doc/qt-4.8/designer-using-a-ui-file.html
        # #widgets-and-dialogs-with-auto-connect
        settrace()
        self.setupUi(self)

How do I fix this?

Windows restart does not help. Deleting .qgis2 folder does not help either.


texts: Google Maps API setTilt(45) on custom tile layer

texts:
With google maps api, I'm trying to put a custom XYZ service on top of google satellite 3d view (45°). Is it possible?

I have a mapnik stack server serving a custom tile XYZ service (transparent polygons) to be consumed by a google maps api application using:

this.tileLayer = new google.maps.ImageMapType({
  getTileUrl(coord, zoom) {
    return [this.tileURL, zoom, '/', coord.x, '/', coord.y, '.png'].join('');
  },
  tileSize: new google.maps.Size(256, 256),
});

This service is working for scales from 9 to 18 and google maps 3d view is set to work at scale 18. But when I reach scale 18, gmaps api automatically switches to 3d view (setTilt(45)) and starts requesting the custom tiles at scale 18, but they're not displayed. Although, google maps satellite images are displayed well and in 3d.

I don't know why but gmaps api requests non-existent tiles, basically having wrong Y value. For example, it requests "18/140819/99179" instead of "18/140819/85969". Therefore, there is not tile to present.

Is there any way to program gmaps api to do on-the-fly tilting (45°) and display a custom tile service?

Ps.: I tried to force the getTileUrl to return this specific fixed tile, but gmaps api doesn't tilt it.

getTileUrl(coord, zoom) {
  if (zoom === 18 && coord.x === 140819 && coord.y === 99179) {
    return this.tileURL + '/18/140819/85969.png';
  }
  return [this.tileURL, zoom, '/', coord.x, '/', coord.y, '.png'].join('');
}


texts:
I'm using Mapbox's Show and Hide Layer example where the user can toggle on/off layers, but seem to be running into an issue I can't pinpoint.

What I did

In Mapbox's example, I replace the first layer (Museum's) with my data and left their "contour" data alone.  My data is an obvious big black blog for simplicity.

The result

The toggle switch appears to work by changing color when turning on/off,  however my big black blog on my map isn't disappearing.

Actions taken

I tried countless times to slowly piece together portions of the code to see where it fails, but just not understanding where I'm going wrong. My guess it's the map.addSource versus map.addLayer which I may not be understanding the difference.  My data was uploaded to a Tileset via a zipped shapefile.
 Perhaps the example being used is from a GeoJson?

Here's my code, except the access token is missing for security purposes.