texts:
So I am currently working on a QGIS Python plugin and I want to update a QDialog with textoutput while an import of data from one file the another is running.

So while it is grabbing data from one file and inserting it into the other I want my output in the dialog too look like:

Importing WMS: connection-wms/example/url: abc.com...
Importing WMS: connection-wms/example2/url: abc.com...
Importing WFS: connection-wms/example3/url: abc.com...

and so on...

My problem is that when I call exec_() on the dialog everything else freezes, like it is waiting for the dialog to do something. When I close the dialog the import runs as intended.

I have allready tried QCoreApplication.processEvents() but somehow can not get it to work.

My current attempt is to run the dialog in its own thread, but with Thread from the threading lib this does not seem to work.

Here is the logic of the import action:

class DatasourceDistributor:

def importSources(self):
    self.dialog = ProgressLog()
    self.dialog.exec_()

    self.dictionaryOfCheckedSources = {**self.dictionaryOfCheckedDataBaseSources,
                                       **self.dictionaryOfCheckedWebSources}

    if self.dictionaryOfCheckedSources:
        self.parser.clear()
        self.parser.read(self.sourceQgisIniFile)

        targetParser = configparser.ConfigParser()
        targetParser.optionxform = str
        targetParser.read(self.targetQgisIniFile)

        for key in self.dictionaryOfCheckedSources:
            iterator = 0

            if key in self.availableWebSources:
                if self.parser.has_section('qgis'):
                    for element in range(len(self.dictionaryOfCheckedWebSources[key])):
                        self.importWebSources(iterator, key, targetParser)
                        iterator += 1
            else:
                for element in range(len(self.dictionaryOfCheckedDataBaseSources[key])):
                    self.importDBSources(iterator, key, targetParser)
                    iterator += 1

        with open(self.targetQgisIniFile, 'w') as qgisconf:
            targetParser.write(qgisconf)

    self.profileManager.updateDataSources()

I feed the dialog with text in another for loop in self.importWebSources(iterator, key, targetParser)

And here the dialog class:

class ProgressLog(QDialog):

def __init__(self, *args, **kwargs):
    super(ProgressLog, self).__init__(*args, **kwargs)

    self.setWindowTitle("Importing Datasources...")
    self.setFixedHeight(300)
    self.setFixedWidth(500)
    self.layout = QVBoxLayout()
    self.progesslog = QTextBrowser()
    self.layout.addWidget(self.progesslog)
    self.setLayout(self.layout)

Please keep in mind that this is my first contact with Qgis Plugins and also with Python

Thanks!


texts: Does anybody know how to only label features of a shapefile that are inside of another shapefile's features used as coverage in Atlas?


texts: can Geocat Bridge enable me to connect to existing PostGIS database and make edits through ArcGIS?

texts:
We are planning to set up a Webgis to share our GIS data on the web. Currently, we are on a desktop GIS and using ESRI's File geodatabase format. We are considering ArcGIS Server or using Opengeo suite. We want to retain using our ArcGIS Desktop but we plan to move our data to a Postgres/PostGIS database to enable multi-user viewing and editing when served on the web.

To be able to connect to a Postgres/PostGIS database without ArcSDE license, there are a number of solutions available, ST-links, GISquirrel, Geocat Bridge...

Geocat Bridge seems appealing cause we can directly publish mxd to Geoserver. But, it's not mentioned on their site whether we can connect to existing  PostGIS data. There is a mention that we can load TO PostGIS but not a thing mentioned on whether we can load data on ArcGIS Desktop from  PostGIS. Anyone who has tried Geocat Bridge? Is this possible in Geocat Bridge?


texts: Creating stacked Python label expression in ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.3.

I have a number of polygons representing exploration licences and a series of attributes including licence name, operator, partners and associated interests (%) within each licence. I am trying to create a stacked label expression that will display in the following order: Licence name, operator and interest, partner(s) name and interest. The way the attribute table is setup is that the percentage interest values are 0 where there is no partner. The dataset is read only and so I can't edit the table to remove the zeros. So far I have the following code:

def FindLabel ([AREABLK] , [Own1] , [Asgnpct1] , [Own2] , [Asgnpct2] , [Own3] , [Asgnpct3] , [Own4] , [Asgnpct4] , [Own5] , [Asgnpct5]):
    if [Own5] is not None:
        return [AREABLK] + "\n" + [Own1] + " (" + [Asgnpct1] + "%" + ")" + "\n" + [Own2] + " (" + [Asgnpct2] + "%" + ")" + "\n" + [Own3] + " (" + [Asgnpct3] + "%" + ")" + "\n" + [Own4] + " (" + [Asgnpct4] + "%" + ")" + "\n" + [Own5] + " (" + [Asgnpct5] + "%" + ")"
    elif ([Own4] is not None and [Own5] is None):
        return [AREABLK] + "\n" + [Own1] + " (" + [Asgnpct1] + "%" + ")" + "\n" + [Own2] + " (" + [Asgnpct2] + "%" + ")" + "\n" + [Own3] + " (" + [Asgnpct3] + "%" + ")" + "\n" + [Own4] + " (" + [Asgnpct4] + "%" + ")"
    elif ([Own3] is not None and [Own4] is None):
        return [AREABLK] + "\n" + [Own1] + " (" + [Asgnpct1] + "%" + ")" + "\n" + [Own2] + " (" + [Asgnpct2] + "%" + ")" + "\n" + [Own3] + " (" + [Asgnpct3] + "%" + ")"
    elif ([Own2] is not None and [Own3] is None):
        return [AREABLK] + "\n" + [Own1] + " (" + [Asgnpct1] + "%" + ")" + "\n" + [Own2] + " (" + [Asgnpct2] + "%" + ")"
    elif  [Own2] is None:
        return [AREABLK] + "\n" + [Own1] + " (" + [Asgnpct1] + "%" + ")"

This gives the following output:

The first few rows work fine but where the partner isn't present it still shows the 0 value with the (%) text added.

I would therefore like to ignore the 0 values where there is no partner in the expression.


texts:
I was given an .mxd that contained icons I need to save and reuse in a webmap.

Is there any way to save these icons as .png or .jpegs instead of .style files?


texts: USA County level, detailed OpenStreetMap basemap in tiff?

texts:
Is there any way to get detailed (close zoom) OSM tiffs at the county level in the USA?

I would like to use them as basemaps for a class I am teaching in GIS, using QGIS.

I know about (and love) the OpenLayers plugin, but it can be flaky (http://hub.qgis.org/issues/5827) and it depends on a decent network connection.

I don't care if the maps are a year or so old, and I also don't care if they are really big.  I would like full coverage of Thurston County WA with full detail, though.

I would most like to have an FTP site to download from, not some app or on-the-fly generated thing.

EDIT:

Here is the workflow I am using, which I think will satisfy my needs:

Download pbf file for my state from Geofabrik
Use Maperitive to generate an MBtile for the region I want (a few square miles only currently, for demonstration, will try the county later today).  I can run this overnight if necessary (generat-mbtile minzoom=1 maxzoom=18)
Load the mbtile into QGIS as a raster format (yay GDAL!)

Thanks to everyone for your help, even though I didn't take any of your ideas in total they helped me figure out the pieces.

I keep thinking I want to produce an "atlas" (metaphorical) of US Census tracts, and one of the things I would include for each tract would be an mbtile and tiff, I think.

It might be a good service to provide these, maybe at the Census Tract level too, for us desktop mappers...

If I am missing something critical, please clue me in....


texts:
Sorry if I sound like I don't understand, but I've been reading around trying to figure out how to do geospatial plotting in Python. I keep running into errors despite the examples I've seen appearing rather within reach.

Currently I have some San Francisco data, where each row is an observation that also contains long/lat data. I would like to do a heatmap of the frequency of the observations. I've come across people using various packages, and I believe some things that are required are shapefiles.

However, for some reason I can't seem to decipher the precise order from which I do things. Does it work like so:

Find a shapefile for San Francisco (with neighborhood polygons)
Plot the shapefile
Overlay a heatmap of the long/lat data

Is this correct? I'm looking for some guidance for a simple way where I can do this using Python 3.x. I've tried packages like Fiona, basemap, GeoPandas, etc., but I'm getting frustrated by not being able to find a simple solution that works in my notebook.


texts:
I have a problem when I start QGIS.

I installed version 3.10 and when I tried to launch it I got the first error window (the one on the left).
I tried to install previous versions (3.8, 3.4) but nothing has changed.

Currently, if I try to launch the program (version 3.8) the window on the right appears first and then the window on the left.

What can I do?


texts: Since ArcGIS 10 comes with the ArcPy package, I'm wondering it is possible to use ArcPy functions to get Symbology (i.e color, width...) of a layer?


texts: Using Creative Commons licenses for spatial data?

texts:
As more people and organizations are realizing the importance and benefits of publishing their data under an 'open' license, people are looking for the best license to apply.

Creative Commons CC has provided licenses that have facilitated great changes in the way that people publish their writing and artwork.

I know that there is some controversy over the use of CC licenses for spatial data, what are the technical concerns?


texts: eCognition segmentation: How to calculate relative area in super-object and relative position to super-object?

texts:
In eCognition, I'm attempting to assess my segmentation results against a set of reference image objects that I have manually digitized.
After performing an overlay function between the segmented image objects and the reference image objects, I would like to calculate the metrics relative area in super-object and relative position to super-object.
I have read in published papers that these metrics can be implemented in eCognition, but I can't work out how to calculate these metrics using eCognition.
Is there a particular tool or algorithm I need to use?
Many thanks,
Sarah


texts: Calculating Area solar radiation gives ERROR 000732: output raster considered as Input?

texts:
I am trying to calculate Area solar radiation for rasters in a list.

It seems that the tool allows to create unique rasters for direct, diffuse and direct radiation duration, and calculates outGlobalRadiation as the final output.

But, when I define the names for the direct, diffuse and direct radiation duration rasters, I got back an error RuntimeError: ERROR 000732: Input Raster: Dataset dir_1 does not exist or is not supported meaning, that my optional output file is actually considered as Input Raster, and does not exist.

Do you understand this error? Or is there any logical error in my script?

import arcpy, os
from arcpy.sa import *

arcpy.CheckOutExtension("spatial")

# Get path
inWD = "C:/Users/2019_bioESS_coord"

# Set working environment
arcpy.env.workspace = os.path.join(inWD, "output/bufRastTwins.gdb")
arcpy.env.overwriteOutput = True

# List all rasters
rasters= arcpy.ListRasters("r_*")

# Calculate solar radiation with defined optional outputs
for raster in rasters:

    # Define names for output optional rasters
    outNameDirect = "dir_" + raster.replace("r_", "")
    outNameDiff = "diff_" + raster.replace("r_", "")
    outNameDirDur = "dirDur_" + raster.replace("r_", "")

    # Define parameters for solar radiation
    in_surface_raster = raster
    latitude = ""
    skySize = 200
    timeConfig = TimeMultipleDays(2018, 100, 150)
    dayInterval = 10    # default 14
    hourInterval = 12   # default 0.5
    each_interval = "NOINTERVAL"  # "INTERVAL"
    zFactor = ""
    calcDirections = 32
    zenithDivisions = 8
    azimuthDivisions = 8
    diffuseProp = 0.3
    transmittivity = 0.5

    # Set names for optional output rasters, seems not working???
    outDirectRad = arcpy.Raster(outNameDirect)
    outDiffuseRad = arcpy.Raster(outNameDiff)
    outDirectDur =  arcpy.Raster(outNameDirDur)

    # Execute AreaSolarRadiation
    outGlobalRad = arcpy.sa.AreaSolarRadiation(in_surface_raster, latitude, skySize, timeConfig,
                                               dayInterval, hourInterval, each_interval, zFactor, "FROM_DEM" ,
                                               calcDirections, zenithDivisions, azimuthDivisions, "UNIFORM_SKY",
                                               diffuseProp, transmittivity,
                                               outDirectRad, outDiffuseRad, outDirectDur)

    print("area radiation saving")
    # Save the output
    outGlobalRad.save(raster.replace("r_", "globRad_"))


texts:
I am trying to figure out the best solution to hide overlapping lines.  I have a layer of polygons representing crustal plates.  I have a polyline layer of fault lines that overlap some of the polygon boundaries.  The fault lines are represented by dotted lines and the polygon boundaries which have solid lines are showing through the dotted lines.

How can I make the boundary lines not show through the polylines?

Currently I can either turn off the borders on the polygons which causes additional problems or I can offset the polylines which is not an accurate representation of the line.

There doesn't seem to be a solution that satisfies my client entirely.  Basically the image shows a similar situation.  The two blue polygons are bordered in white and the fault line in dotted black which he doesn't want to have white behind.  Basically he wants to have no border where there is a dotted line which is impossible.  My advise to him is to either remove the borders all together or live with the overlap.


texts: Should I use geometry type over geography type when usually expecting back a lot of rows?

texts:
I am in the planning stages of a service that will do a fair share of searches for items in a geographical database.
I will have to use PostgreSQL and PostGIS for that purpose and have learned so far from this article that geography searches returning a lot of rows come with a high performance penalty.
My service will be based around searches asking for every item within a radius from 1km up to 50km.
So I will have to utilize the ST_DWithin function a lot.
I am expecting to get back several hundred rows per square km in urban areas but in rural areas probably only a handful.
[edit: added after discussion in the comments]
This made me think so far I might be better of building everything right away in geometry types, until I read in this answer to a question about ST_DWithin() that the use_spheroid parameter=false was added to lower some performance issues. But I could not find out if this would be enough to decide to use geography.


texts: FME: Identify a string within a field and populate another field with that same string

texts:
I have a shp which needs to be split according to attributes in a field. The problem is that the field contains other Characters and dates. For Example the Field is called "Landtype" and contains the following "Woodland_20120502". The Date changes for each data record so I can´t just select all records with this value. I need to first select all records that have the word "Woodland" and then populate another field with the word "Woodland" in FME or ArcGIS modelbuilder or Python.

Could anyone advise me on how I would best do this?


texts: GeoPandas returning incorrect spatial join results. GDAL problem? GEOS problem?

texts:
I'm having a problem that I can reproduce on an Ubuntu installation of GeoPandas (0.3.0), but not on my local OS X intallation (same geopandas version). Here's the example (taken from the geopandas documentation):

import geopandas as gpd
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))
countries = world[['geometry', 'name']]
countries = countries.rename(columns={'name':'country'})
cities_with_country = gpd.sjoin(cities, countries, how="inner",  op='intersects')
cities_with_country.head()

The correct result, as shown here and as obtained on my OS X installation, is:

             name                                     geometry  index_right  \
0    Vatican City  POINT (12.45338654497177 41.90328217996012)           79
1      San Marino    POINT (12.44177015780014 43.936095834768)           79
192          Rome    POINT (12.481312562874 41.89790148509894)           79
2           Vaduz  POINT (9.516669472907267 47.13372377429357)            9
184        Vienna  POINT (16.36469309674374 48.20196113681686)            9

     country
0      Italy
1      Italy
192    Italy
2    Austria
184  Austria

However, on my Ubuntu install, I get this:

      name                                     geometry  index_right  \
183  Kabul  POINT (69.18131419070505 34.51863614490031)            0
183  Kabul  POINT (69.18131419070505 34.51863614490031)            0
183  Kabul  POINT (69.18131419070505 34.51863614490031)            0

         country
183  Afghanistan
183  Afghanistan
183  Afghanistan

and len(cities_with_country) is only 3. The same problem happens with other data as well.

I've check all the listed dependencies on both installs and verified that they're the same versions. However, the GDAL versions are different. The correctly functioning OS X install has:

$> pip2 freeze | grep GDAL
GDAL==1.11.5

While the incorrectly functioning Ubuntu install has:

$> pip2 freeze | grep GDAL
GDAL==2.2.2

I'm reluctant to start messing with downgrading my GDAL installation, so I'm hoping that someone can confirm or refute my theory that geopandas is (indirectly) dependent on GDAL and that version 2.2.2 could be causing the problem. ...and I'm wondering if downgrading is the only solution.

After revisiting this and looking through the geopandas and shapely code, it seems more likely that this would be related to GEOS rather than GDAL, right? As far as I can tell, both installations are using GEOS from the wheel installed with shapely. Both give the same geos version string:

>>> import shapely
>>> shapely.geos.geos_version_string
'3.6.2-CAPI-1.10.2 4d2925d6'


texts: Finding route contains Source and Destination points

texts:
I have a problem with finding routes that are near (in range X km) the passed points (for e.x. A and B) in PostgreSQL PostGIS.

Example:

The point A denotes Source point and The point B denotes the Destination.

Tn system there are routes represented as LineString. I am also able to consider another ways to represent the routes.

Phere is no problem to find routes contains the points A and B in specific range (for e.x. 5 km):

SELECT * FROM routes
WHERE
    ST_Distance_Spheroid(route.lineString, 'POINT (long_A lat_A)', 'SPHEROID["GRS 1980",6378137,298.257222101]') < 5000
AND
    ST_Distance_Spheroid(route.lineString, 'POINT (long_B lat_B)', 'SPHEROID["GRS 1980",6378137,298.257222101]') < 5000

But there is the problem, that this query considers the path regardless its direction. It is insensitive to order the points that describes the linestring.

So there is a question. How select the routes contains points in specific range and consider order of points.

BTW: There is possible to represent routes as points and store in database the POINT and order number. So I can look for nearest points.


texts: About support vector machine classification

texts:
I want to classify Sentinel 2 images with the SVM in orfeo toolbox in QGIS.
I have made a ROI (2 classes) and use TrainImages classifer (libsvm). I have made a SVM model,but when I use ImageClassification as a results I get a gray images.
Can you help me about this problem. I try evrithing but always I get a gray images.It is a images with a shades of gray. I thinkg that OTB do classification but I dont know why images is gray


texts: PgRouting: What index should have the restriction table

texts:
I have this tables to feed my routing engine using pgr_trsp:
vzla_rto have the links and nodes (650k)
vzla_restriction have the turn restrictions (10k)

CREATE TABLE map.vzla_restrictions
(
  rid integer NOT NULL,
  to_cost numeric,
  to_edge integer,
  from_edge integer,
  via text
);
CREATE INDEX vzla_restrictions_from_edge_idx ON map.vzla_restrictions USING btree (from_edge);
CREATE INDEX vzla_restrictions_to_edge_idx ON map.vzla_restrictions USING btree (to_edge);

I was checking table statistics and realize i put an effort to create index for vzla_rto, but looks like i didn't enough for vzla_restrictions.

                    | seq_scan      | seq_tu_read     | idx_scan        | idx_tup_fetch | n_tup_ins
vzla_restrictions   | 1,672,959     | 16,848,370,089  |      -          | -             | 10071
vzla_rto            | 1,549         | 977,994,776     |  862,610,399    | 3,866,967,927 | 648968

So which index should i add to vzla_restrictions to help pgr_trsp()?


texts: Using Set Z-Values in QGIS Processing Modeler

texts:
I have a set of contours that I am trying to automate a process on. This process clips the contours by a boundary, sets z-values based on the elevation field, and then exports those contours as a dxf file to be used in CAD. My issue is that I cannot figure out how to set the Z-value from the Processing Modeler algorithm.
Is there a certain syntax that I need to use so that I can assign the "ELEV" field from the output of my clip function to the Z-value in the Set Z-value function? I see no drop down or option to select a field from the output of the clip function so I am confused as to how to accomplish this.


texts:
I am using the intersect service from the ArcGIS Geometry service.
referring : http://help.arcgis.com/en/webapi/sil....htm#Intersect

I have modified it to make it work with polylines
it works fine to intersect a polyline with a polygon but does not return anything when I try intersecting a polyline with a polyline?


texts: Using QGIS to convert from WGS to UTM not working

texts:
My predicament is this. I have inherited a bunch of data in multiple projections, and I need to calculate hub distance to get distance from one grid of points to a contour layer (I have extracted vertices for this layer so it is just distance points-points).

The first layer when looking in the metadata is in WGS84 - pretty standard. The other grid of points is in UTM (GDA 94, 55 G).

The confusion is this:

When I open both files in a single workspace with no specified projection (the projection in lower right is listed as NO PROJECTION, with a bunch of extra code afterwards), I can see both files and they are more or less arranged how I would have assumed. However, this means when I calculate hub distance I get degrees - I need kilometres.

SO, I try to save the contour file that is in WGS84 in GDA94. This works fine, but if I open a new workspace and set project CRS to GDA94 and open the grids, I see them where expected at the coordinates I expects, however when opening the contour file it becomes flipped and distorted, and placed thousands of kilometres from where is was when I had them open in default CRS - presumably projected on the fly. My question is this - if they both look normal before I export them... What am I doing wrong for the contour points to be so drastically warped and moved??


texts: How to rename large amounts of spatial data?

texts:
I have approximately 1500 raster files following this format:

m_3609922_nw_14_1_20100425

I need these subsetted to, for example: 3609922nw.  What are appropriate methods to rename large numbers of spatial data, that may or may not contain ancillary files (e.g. .xml, .ovr)?


texts: Spatial functions with MySQL, XAMPP, and Wordpress

texts:
I am trying to figure out how to using the ST_AsGeoJSON function in MySQL, but I keep getting an error that says that this function does not exist. I have found in the MySQL documentation that this spatial function does exist.

mysql> SELECT ST_AsGeoJSON(ST_GeomFromText('POINT(11.11111 12.22222)'),2);

https://dev.mysql.com/doc/refman/8.0/en/spatial-geojson-functions.html

My current setup is Windows 10 with XAMPP with Wordpress.


texts: Improving Geostatistical Analyst Wizard performance?

texts:
I have performed CoKriging using Geostatistical Analyst in ArcGIS 10. When I tried to export the prediction to raster, the process became so slow.

Is there a way to make the exporting faster?


texts:
I currently have a line that I want to but a label onto in Tilemill. It currently looks like this:

But what I would prefer is something like this:

where the line isn't rendered around (and especially underneath) the label. Is there any way to archive that in Tilemill?


texts: CartoDB Runtime Visualization on Local Instance

texts:
I'm trying to create a runtime visualization on a local instance of CartoDB. Here's the relevant code:

cartodb.createLayer(map, {
            user_name: 'testuser2',
            type:'cartodb',
            tiler_host:'', //??
            tiler_port: '', //??
            tiler_protocol: '', //??
            sublayers: [{
                sql: "SELECT * FROM metro_stop AS pop",
                cartocss: "#pop {marker-fill: #0fb16a; marker-width: 10;}"
            }]
        }, {no_cdn: true}).addTo(map);

The problem is, I don't know the values of the three missing parameters: tiler_host, tiler_port, and tiler_protocol. Does anyone know how I would find those out?

Edit to add:
Here's the error I'm receiving when I set tiler_host to 'localhost.lan', tiler_port to '8181', and tiler_protocol to 'http':

http://testuser2.cartodb.com:8181/api/v1/map?config=%7B%22version%22%3A%221.0.0%22%2C%22stat_tag%22%3A%22API%22%2C%22layers%22%3A%5B%7B%22type%22%3A%22cartodb%22%2C%22options%22%3A%7B%22sql%22%3A%22SELECT%20*%20FROM%20metro_stop%20AS%20pop%22%2C%22cartocss%22%3A%22%23pop%20%7Bmarker-fill%3A%20%230fb16a%3B%20marker-width%3A%2010%3B%7D%22%2C%22cartocss_version%22%3A%222.1.0%22%7D%7D%5D%7D&callback=jQuery17209372722392436117_1401785836614&_=1401785836710
Failed to load resource: Could not connect to the server.

For whatever reason, this is pointing to CartoDB.com to try to get the info.


texts: Saving a WFS layer directly to the disk using QGIS or PyQGIS?

texts:
I know how to properly load a layer in QGIS from a WFS service, but I would like to directly save it to a .shp file (without the preliminary loading of it into QGIS).

Is it possible to do this within QGIS or with a PyQGIS script?

A sample URL from which starting would be this one:

http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/wfs/Carta_geologica.map

having these capabilities.


texts: QGIS log Error displayed within Quick OSM Plugin

texts:
I am having an issue using the Quick OSM Plugin. I am using QGIS V3.4. I am following this tutorial linked below.

QGIS Training Manual

The issue I am facing is occuring at Step 9 of section 1.2.1 'Try Youreself Create OSM based Vector Files. After downloading the sample data used within the 'QGIS Traning Manual'

I get this error message --> "Error in the QGIS Logs, QuickOSM Panel, please report it to GitHub". I have reported this issue, however was wondering if this is my human error or if the software or data has the issue and can any one solve this?


texts: Turn off confict resolution in GeoServer 2.6.2

texts:
Labels are being hidden in GeoServer on a single tile image. I suspect this is because they are too close to one another.

Based on other questions here, I have tried to add:

<VendorOption name="conflictResolution">false</VendorOption>

to the TextSymbolizer. It doesn't seem to do anything. Neither does spaceAround with a negative value. At the moment I am not displaying the polygon, just the labels (although I would like both).

I want to turn conflictResolution off completely, and show all labels, even if they overlap.

These two parameters are mentioned in this question, where they appear to have solved the persons problem: GeoServer VendorOption for SLD to place labels overlapping and out of bounds

(they have another problem of labels being near the edge of the image, which I do not have).

Here is my SLD:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <!-- a Named Layer is the basic building block of an SLD document -->
  <NamedLayer>
    <Name>area</Name>
    <UserStyle>
      <FeatureTypeStyle>
        <Rule>
          <TextSymbolizer>

           <Label>
             <ogc:PropertyName>shortest_label</ogc:PropertyName>
           </Label>
            <Font>
                <CssParameter name="font-family">Arial, Verdana, Sans-Serif</CssParameter>
                <CssParameter name="font-size">12</CssParameter>
                <CssParameter name="font-weight">normal</CssParameter>
            </Font>
            <Halo>
              <Fill>
                <CssParameter name="fill">#ffFFFF</CssParameter>
              </Fill>
            </Halo>
            <VendorOption name="conflictResolution">false</VendorOption>
            <VendorOption name="spaceAround">-1</VendorOption>
          </TextSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts: Calculating spatially weighted average value of vector layer in QGIS 3

texts:
I have a vector layer across an area with several different polygons of varying size. Each polygon has a rate of fertilizer as an attribute. I would like to calculate the average rate of fertilizer across the area.

Initially, I used the tool "Basic statistics for fields" but I realized that this tool just takes the average of the numbers in the attribute table without adjusting for the polygon sizes. This gives an incorrect average rate since there could be a large polygon with a very low rate, however, this large polygon would count towards the average as much as a very small polygon.

Is there a tool that allows me to calculate the average rate of fertilizer so that when I multiply the average rate by the area, I would get the correct amount of fertilizer needed for the farm?

I can calculate the correct average in Excel or using the field calculator, but I think there must be a faster way to do this in QGIS 3.12.0.


texts: Determining off to add scale bar and north arrow in spplot

texts:
I have a spplot  in which I want to add north arrow and scale bar.
bbox(dis) # dis is my shapefile, an sp object

        min     max
x -285706.6 1492871
y 2620690.6 4120742

I want to attach the north arrow and scale bar. My project UTM 42N


texts:
What does LW stand for in liblwgeom of PostGIS?

I could not find any explanation anywhere, or in the source code (mirrored) at https://github.com/postgis/postgis/tree/svn-trunk/liblwgeom


texts: Finding clusters of points based distance rule using R?

texts:
I have some points that I need to divide into groups based on clusters of points:

x <- c(-1.482156, -1.482318, -1.482129, -1.482880, -1.485735, -1.485770, -1.485913, -1.484275, -1.485866)
y <- c(54.90083, 54.90078, 54.90077, 54.90011, 54.89936, 54.89935, 54.89935, 54.89879, 54.89902)

library(sp)
xy <- SpatialPoints(matrix(c(x,y), ncol=2))
proj4string(xy) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
xyTransformed <- spTransform(xy, CRS("+init=epsg:27700 +datum=WGS84"))

plot(xyTransformed)

Each group should be defined using this rule: all points in the group should be within 40m of each other.

Based on the points in code above, I believe I should end up with 3-4 groups.

Does anybody know how to do this in R?


texts: Why does no data show up in the embedded CartoDB map?

texts:
It seems something went wrong with CartoDB. Recently I am trying to embed my visualization into my website. But only comes with a blank map. Then I noticed that the live example from the official website also has the same issue.

Here is the link:
http://cartodb.github.io/cartodb.js/examples/easy.html


texts: Setting legend area proportion on webmap produced by qgis2web for a mobile phone/tablet

texts:
I have produced a webmap using qgis2web tool using the open layers option. It is hosted on a remote server. Viewing it on a desktop/laptop is fine. However when i try to view the same using a phone,i find that the legend is taking up most of the screen space!

How do i edit the code to allow the legend to take up specific percentage of available screen space?

@image is how it is displaying on a mobile phone


texts: Rasterizing Shapefiles onto TIFF loses style

texts:
I have two contour shapefiles, one for 2ft and one for 10ft contours.  I want to put these on a TIFF file I have of the aerial.  I achieved this by going to Raster > Conversion > Rasterize and using the shapefile as input and the TIFF as the output and when this happens I get contours of a red and turquoise color.  I want the contours to be the same color all throughout the TIFF I load them on and also I want to be able to choose the color and style that I use.

I have tried changing styles of the contours and re-saving and also changing the Band type of my raster since it is a multi-band and I figured this was the source of my problem.  However, whenever I try this the TIFF file becomes different colors and unusable.  Attached is a picture with two arrows pointing to the two different colors the contours come in when I go through the Rasterize tool.


texts: QGIS distance between centroids from NUTS1 information

texts:
I have calculated the centroids from NUTS1 information. Now for every neighbor centroids, I would like to calculate the distance. I want to connect only the neighbor's Centroids regarding the NUTS1. I tried the formula join attributes by nearest but as the title suggests joins the nearest and not the neighbors.
Here is a screenshot of the problem

Any suggestions?
Furthermore, after tries, I manage to solve the problem. First I created a distance matrix with k nearest equal to 7. Then I duplicate the distance matrix and by using "Join by lines" I connected all the relevant centroids. At this point, I start manually deleting the lines that were not correct according to the neighbor's idea. After that, I checked the attribute table (in order to save in a csv) but there were multiple duplicates for the same connected line. Perhaps the command Join by lines caused the duplicates problem. I export it in excel. Where I could erase easily all the duplicates with a command.
The map looks like that right now


texts:
I have two sets of data.  One comprises all home sales in a given area, where each data 'point' is actually a polygon representing the lot that was sold.  The other set comprises arrests made in the same area, where each point also has a date associated with it.

I want to use the IDW tool to tell me roughly how much arrest activity occurred around a given home in the 12 months preceding the sale.  I essentially want to go row by row in the house sale data and tell the spatial analyst tool to consider only points from the second dataset from date n to n-365, then interpolate all points within 1000 feet that satisfy this condition.  Any ideas?

I believe I can select all the points that meet my date criteria, but I am not sure how to feed that to the IDW tool interactively (iteratively - i.e. a different selection criteria for each data row, because each home sale happened at a unique location on a unique date).


texts: Difference tool QGIS shows too many Differences

texts:
After running QGIS Difference Tool, I get a lot of differences, but I don't know where they are coming from. The input files are almost similar...

Input Image (overlay input almost the same, but from different years with one or two different polygons):

Difference Output:


texts: Using Modelbuilder or Python to write output a folder name using database or table file?

texts:
Has anyone worked on Modelbuilder or python to write to automate creating a folder using the database or table file to write output a folder name and or/create the subfolders? I recalled I had an aml that allowed me to modify codes to allow to output new folders within the folder.

For example, you have a database that has the list of the name of the counties on the dbf file or table and then insert it to the Modelbuilder. Example name of county : Dona Ana and then use that name to create a folder Dona Ana  ( to hold imagery for NAIP on the server when I get all of the New Mexico counties in a hard drive to me)

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

I'm adding this screenshot to show what I am doing wrong... ( This is the pre-screenshot)
that was done yesterday December 9th.

Now that that i have figured it out and here is the screenshot December 10th.


texts: Is there a TOC control available for Esri's Silverlight/WPF SDK?

texts: I've got a requirement calling for a TOC control on a Silverlight/WPF application.  Does anyone know of a TOC control (analogous to the TOCControl in ArcGIS Engine) that I can hook into Esri's Map control?


texts: I have a simple interaction on my map, by this interaction user can draw circle in any part of map. But i need to consider a current projection, can i apply projection to interaction? I try to say what if we draw circle on map center is looks like a circle, but if we draw it closed to poles - its looks like ellipse? May be somebody can help me.


texts: Creating layers using definition query for each field in ArcPy?

texts:
I have a shapefile that contains a text field "Company" that i would like to use to run a definition query to create a layer file. I want to be able to create a layer for each company.

I had a look around and I found Creating new layers using ArcPy to define Definition Query for each field value? which seems to detail the process.

In my code i am trying to write out a layer file for each unique value that determines a definition query. My source file contains a shapefile of shops , within the shp file there is a field that contains the company that owns the shop. I want to run a definition query on each unique company contained within the company field

i have been able to get the unique companies into a list( as pointed out , i needed to use a set to get the unique values of the search curosr) , that i would like to loop through. My code is as follows

import arcpy
from arcpy import env
# This is where the file is being picked up from
arcpy.env.workspace = "S:\\AppsData\\GIS\\PROJECT\\COMPANIES\\SHAPEFILES\\"
# The file to be processed
Source_File = "SDE_COMPANY_LIST.shp"

Company_Bulk_Company_List = [row[0] for row in arcpy.da.SearchCursor(Source_File,"COMPANY")]
UniqueCompany = set(Company_Bulk_Company_List)
print UniqueCompany

# Getting the list out

for company_query in UniqueCompany:

env.workspace = "S:\\AppsData\\GIS\\PROJECT\\COMPANIES\\SHAPEFILES\\"

out_layer0 = str(company_query)+"lyr"
in_layer = out_layer0
out_layer = in_layer +"lyr"

#MakeFeatureLayer variables
in_features = "SDE_COMPANY_LIST.shp"

where_clause = '"COMPANY" = \'company_query\''
print where_clause
workspace = "S:\\AppsData\\GIS\\PROJECT\\COMPANIES\\SHAPEFILES\\"

try:
    # Execute MakeFeatureLayer
    arcpy.MakeFeatureLayer_management(in_features, out_layer0, where_clause, workspace)

    # Execute SaveToLayerFile
    arcpy.SaveToLayerFile_management(in_layer, out_layer, "ABSOLUTE")
except:
    print arcpy.GetMessages()

I the layers are being created , but there is a problem with the definition query as what i am in getting inside the defination box for the layer is

company = ' company_query'

where is i like to see
company = ' variable that is inside company_query'


texts: Monteverdi - Orfeo Toolbox: OTB Applications browser empty

texts:
I am on Ubuntu 16.04 and installed Orfeo Toolbox using the binaries[¹]

I loaded Monteverdi and clicked on View -> OTB Applications browser but the Menus is empty:

[¹] https://www.orfeo-toolbox.org/download/

How can I solve this?


texts:
I have a layer (streets) where I want to remove duplicate labels and then set a display interval for the labels. This is relatively easy in ArcMap because removing duplicates is based on the label text, not the geometry. I set the remove duplicates radius to the same same label repeat distance (although in most cases I could get away with omitting the repeat distance) and I'm good to go. In QGIS it's not so easy.

I am using QGIS3 and labeling major highways for my region. Assume a roughly rectangular region slightly higher than wide, with major north-south and east-west arteries bisecting at roughly the centroid of the region. I have tried every combination of settings I can come up with and the best I get is the label for the east-west artery displaying twice west of center (with less than 4" between the labels) and once east of center, and the north-west artery is labeled twice south of center (with less than 4" between the labels) and there are no labels at all north of center. (My Composer layout is 36" wide x 42" high, with the region filling most of that area).

My setting for label Placement is to set Repeat to 4" and for label Rendering I have checked Merge connected lines to avoid duplicate labels.

I found this solution (Removing duplicate point labels in QGIS) but it is for an earlier version of QGIS which uses Python 2.7. If possible I would prefer a solution through the interface but I'm not averse to using a little Python, but not being conversant with QGIS3's Python object model I'd be spending an inordinate amount of time with this; I've already spent a lot of time trying to make these labels work.


texts: Incorrect displaying of icons in ModelBuilder?

texts:
I'm using ArcGIS 10.1 for Desktop.

Something weird happened with icons in my model.

They are displaying incorrectly as you can see here

.

In other models all is just fine.


texts: What formula is used for distance between 2 points in Google Maps?

texts: What formula is used by Google Maps for calculating the shortest distance between 2 points (having their Lat, Lng coordinates)? Is it purely based on Haversine? Or is it something different? Is there a mathematical or code implemenatation reference available somewhere?


texts: If then statement in model builder for several features after a clip

texts: I'm merging several databases and clipping out certain areas of interest.  Sometimes there are no points that are clipped because there is nothing in the clip.  I have each of the features output to a geodatabase.  Is there a way for it not to export the feature class to the geodatabase if the file is empty?


texts: Programming for a beginner (with GIS experience)

texts:
I am a GIS Specialist with an Earth science/geography background, and about 8 years experience using both ArcGIS and QGIS.

I have written small snippets of Python code to help myself with basic tasks such as editing attribute tables or displaying map labels and symbols in interesting ways.

I have also used ModelBuilder tools in ArcGIS to automate tasks and simple work flows. But I would really like to learn more about computer programming in general.

I am planning to take a 24-hour ESRI classroom course on "Geoprocessing Using Python" next month. But what can I do over the next few weeks to get better at programming, in preparation for the class?

I really want an online course or textbook that is for beginners, that focuses on things like loops, if/then/for statements, etc.  I never seem to learn much at all from those "here is the code, type it in, it will do x" tutorials.

I think I really need something that is a combination of theory and practical knowledge. A focus on Python and/or GIS is of course a bonus.

Any suggestions?


texts: Square polygons that are covered by trees

texts:
I would like to know if there is a way to square polygons?

I would like to calculate the total square meters of houses. I have extracted the buildings as well as tree coverage but I am unable to get a full total square meter of the house due to the tree coverage that covers the building.

I was thinking if there is a way that I can square out the building to try and accommodate the tree coverage?


texts: Is it possible to export the Google Earth Engine basemap (satellite version, not map), instead of something from an image collection?

texts:
I was originally going to use Landsat8 or Sentinel Data, but they don't go as clear as I thought they would. I'm looking for something at around 1m per pixel (Landsat8 is about 30m per pixel), like the Planet SkySat data (https://developers.google.com/earth-engine/datasets/catalog/SKYSAT_GEN-A_PUBLIC_ORTHO_RGB); this does not seem to be available for my country.
When I look at the basemap and set it to Satellite, it zooms in enough for what I need. Is it possible to download this instead? The codes I'm currently using are as follows:
var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1')
                  .filterBounds(ee.Geometry.Point(120.95, 14.55))
                  .filterDate('2018-09-02', '2020-09-02')
                  .sort('CLOUD_COVER')
                  .first();

var trueColor432 = dataset.select(['B4', 'B3', 'B2']);
var trueColor432Vis = {
  min: 0.0,
  max: 20000.0
};

var exp_geometry = ee.Geometry.Rectangle(
[120.9519459437584, 14.558496385034239,
121.02730528824083, 14.645375475871372]);

var trueColor432_viz = trueColor432.visualize({
    min: 0.0,
    max: 20000.0
});

Map.setCenter(120.95, 14.55, 10);
Map.addLayer(trueColor432, trueColor432Vis, 'True Color (432)');

Export.image.toDrive({
 image: trueColor432_viz,
 description: 'Data_Landsat8',
 scale: 30,
 region: exp_geometry,
 folder: 'Data'
});

As seen, they start with the ee.ImageCollection() function, which I'd like to skip.


texts: Convert xy spatial points into spatial lines

texts:
I have produced a SpatialPoints object (named v3) in R from the function gIntersection (package rgeos). The SpatialPoints object v3 has the following structure and is represented by red points on the image below.

    > summary(v3)
    Object of class SpatialPoints
    Coordinates:
            min       max
    x  310440.5  310680.5
    y 4996427.5 4997297.5
    Is projected: TRUE
    proj4string :
    [+proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
    +units=m +no_defs]
    Number of points: 32

From the SpatialPoints object v3, how can I convert the red points into 3 spatial lines as in the image below ?

The SpatialPoints object v3 has only spatial coordinates x and y.

Update

I used the function gIntersection with two objects of class spatialPolygons (squares in the images) and I almost get the desired result (i.e., to draw boundaries between polygons that have different values). But the function  gIntersection creates a object of class SpatialCollections with 1 feature. So, How can I convert objects of class SpatialCollections into objects of class SpatialLines ? From the reproducible example, I would like to get only the lines in red (and not the points).

Here the reproducible example:

r <- raster(ncols=90, nrows=50, resolution = 30)
projection(r) <- CRS("+proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
values(r) <- sample(1:3, ncell(r), replace=TRUE)
plot(r)

p1 <- rasterToPolygons(r, fun=function(x){x==1})
plot(p1,add=T)
p1 <- as(p1,"SpatialPolygons")

p2 <- rasterToPolygons(r, fun=function(x){x==3})
plot(p2,add=T)
p2 <- as(p2,"SpatialPolygons")

p3 <- gIntersection(p1,p2,byid=T)
plot(p3,add=T,col="red",lwd=3)
> p3
class       : SpatialCollections
features    : 1
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0


texts:
I'm using QGIS Desktop 3.10.3 with GRASS 7.8.2

I tried to run the plugin "Data Plotly"
https://anitagraser.com/2017/12/06/data-exploration-with-data-plotly-for-qgis3/

After installing the plugin I get the following message:

I already tried to reinstall the plugin and install also an older QGIS Version but still can't run the plugin.
I am thankful for every help. Unfortunately, my scripting skills are not enough to fix the problem on my own


texts:
I'm in a introduction course for ArcGIS and I'm just starting to work on my project. I'm trying to create a map showing the relationship between poverty and education. I also want to try and present other correlations like race, healthcare, access to services, environment.

My question is regarding the FactFinder website for the Census Bureau. How do I get data for household income for each state? Should I use the FactFinder page? Are there other ways of obtaining this data?


texts:
I'm using arcmap 10.1 and have one shp with many polygons with different typeID. I want to make a serie of maps showing each typeID on an individual map. Since I have many typeIDs I want to use arcpy. And here is the problem. I don't know which arcpy function to use to select only one typeID and remove all the others; and then repeat this step for all IDs.

SOLUTION:
Here is my code after I fixed it. Now it works:)

mxd=arcpy.mapping.MapDocument("CURRENT")

def getValueList (inputTable, field):
    values = set()
    rows = arcpy.SearchCursor(inputTable)
    for row in rows:
            values.add(row.getValue(field))
    return sorted(values)

list=getValueList('my_shp_file', 'typeID')
lyr=arcpy.mapping.Layer('my_shp_file')
where_clause='"typeID"='

for type in list:
    lyr.definitionQuery=where_clause+str(type)
    arcpy.mapping.ExportToJPEG(mxd, "C:\\my_location\\"+str(type)+".jpg")

Thanks for suggestions!
Rok


texts: How to export multiple Shapefiles to one KML file in QGIS?

texts: As the title says, is it possible to export multiple Shapefiles to one single KML in QGIS (Similar to "Map to KML" in ArcMAP)?


texts: QGIS CRS reprojection creates empty files and misbehaves

texts:
The problem mostly centers around reprojection in QGIS 2.4.0 under Windows 8.1. I am trying to reproject several files (world map shp file from NaturalEarth (ne_10m_land), Coastline map (ne_10m_coastline) or csv files) from WGS 84 to World Mercator as I need units in meters for a buffer analysis (global scale).
The world map (originally WGS 84) will reproject when the layer CRS is set to the new CRS (using Set Layer CRS) but trying to do a buffer analysis on the reprojected world map in World Mercator results in empty attribute tables. Same happens when saving the reprojected map in the new CRS (using Save as). The buffer analysis creates a new (empty) file but does not operate at all. Tried several other projected CRSs and it is either the same, the map won't reproject at all or QGIS crashes while trying to do it.

Trying to load the CSV file in a Projected CRS layer results in the layer having the new CRS but points being displayed somewhere completely different, no reprojection to other Projected CRS nor saving in new CRS changes that.


texts: Draw graphic on screen without requiring double click to finish sketch?

texts:
I am trying to get the user to draw a line (only 2 vertices) on the screen (just like the georeferencing link/spatial adjustment link).

I got the below snippet from Esri website which does give me a polyline but it requires the user to double click to finish the sketch.

public IPolyline GetPolylineFromMouseClicks(IActiveView activeView)
        {

          IScreenDisplay screenDisplay = activeView.ScreenDisplay;

          IRubberBand rubberBand = new RubberLineClass();
          IGeometry geometry = rubberBand.TrackNew(screenDisplay, null);

          IPolyline polyline = (IPolyline)geometry;

          return polyline;

        }

I would like the tool to finish the drawing as soon as two points are clicked, thus avoiding a double click.

Can somebody help me on this?


texts: Why is ICalculator.Calculate significantly slower than the CalculateField geoprocessing tool with SDE datasets?

texts:
I just discovered today that there are essentially two ways to calculate field values in ArcObjects:

1. Using ESRI.ArcGIS.GeodatabaseUI.ICalculator:

ITable table = …;
ICursor updateCursor = table.Update(null, false);

ICalculator calculator = new Calculator()
{
    Cursor = updateCursor,
    Field = …,
    Expression = "…"
};
calculator.Calculate();

Marshal.FinalReleaseComObject(updateCursor);

2. Using ESRI.ArcGIS.DataManagementTools.CalculateField:

var calculateField = new CalculateField()
{
    in_table = …,
    field = …,
    expression = "…",
    expression_type = "VB"
};

var gp = new Geoprocessor();
gp.Execute(calculateField, null);

While performance of the two approaches is similar with File GDB tables, I've found that it varies a lot for SDE tables: I've measured a processing speed of 21 milliseconds per row with the former approach (ICalculator), while geoprocessing requires only 10 milliseconds per row (<50 %).

Does anyone know, or have an educated guess, why ICalculator is slower?

With the ICalculator approach, it doesn't seem to matter whether the cursor used is a recycling one or not. But apparently, an update cursor is required. (I tried using a search cursor, too, but that resulted in warnings and errors — discoverable via ICalculatorCallback — immediately after calling Calculate.)
Also with the former approach, I tried reducing network bandwidth usage by specifying only the minimum required SubFields in the IQueryFilter used to retrieve the update cursor. Surprisingly, this doesn't improve runtime performance at all.
I tried all above methods outside any edit session or explicit transaction.


texts: QGIS SPIT errors loading shapefile into PostGIS

texts:
With QGIS 1.8 on Windows 7, I just installed PostgreSQL / PostGIS and tried to load my first data into it, a shapefile using SPIT, but failed:
ERROR: Geometry type (Polygon) does not match column type (MultiPolygon).
Then I used a QGIS Vector tool "MultiPart to SingleParts", which changed the multipolygon type to "polygon" and tried to load that newly generated shapefile, which looks okay in QGIS, using SPIT once again.  This time I got:
ERROR: zero-length delimited identifier at or near """"

My QGIS Help dialog says that I'm using PostgreSQL Client Version 8.3.10, but I had just installed PostgreSQL 9.2 with PostGIS 2.0 using the StackBuilder from Enterprise DB.

I just found that I can load the shapefiles without error using the "PostGIS 2.0 Shapefile and DBF loader exporter" that installed from StackBuilder, and QGIS can then load that layer.

Is there some inherent incompatibility with what I just installed from Stackbuilder (to be able to create PostGIS databases) and is there a simple solution for the QGIS SPIT tool to get it operational?


texts: Add a random value of the order of 10^-2 to raster layer?

texts: I was doing some raster analysis for smoothing in ArcGIS, one of my task is adding a random value of the order of 10^-2 to a raster layer. Does anybody know how to do it?


texts: Displaying reference pictures to a varying value in QField

texts:
While developing a QField project for forestry inventories, I aim at displaying botany board related to the current selected species selected in a scrolling-menu for field operator. It should be used as a reference to compare with the species suspected on the field. I followed this procedure : https://qfield.org/docs/prepare/add-1-n-pictures.html
Here is the referencing layer as a table containing pathes for the pictures related to the species named "photos_essence.csv"

Here is the referenced layer "Qfield_comptage.shp" used for digitizing.

In the project properties, i added a relation "photos_essence" between the referenced layer "Qfield_comptage.shp" and the referencing layer "photos_essence.csv". It should be linking species ("Essence") identified and the path to the pictures related to those species.

I added the relation "photos_essence" in  Attribute Form of the "Qfield_comptage.shp" referenced layer.

Then I changed the widget type to attachment, linking it to the file containing the pictures

The end result gives a null value instead of displaying the pictures related to the species identified.

Do you have any insight of what i'm doing wrong ?
I do not know if such a widget is conceivable.


texts: GRASS raster calculator - fill nodata values from other raster

texts:
I want to use r.mapcalc to calculate a raster from two inputs: water_rst and contour_rst. I need a raster with all values from water_rst and then values from contour_rst where water_rst has no data.

An expression like comb_rst = not(if(water_rst))*contour_rst produces raster with only contour values, but no initial water values. How to modify this expression to get both water values and contours?


texts: Has anyone tried to use the Sentinel toolbox in command-line? I have downloaded the SNAP tool and installed on a Ubuntu workstation and the GUI works well. I now want to want to use it in a command-line so that I can automate some processes. Is there any documentation on how to get started with the SNAP tool on a command-line environment?


texts: Projecting all the points from the csv file at (0,0) in a map

texts:
I am trying to plot some points on a map using the GeoTools tutorial (http://docs.geotools.org/stable/tutorials/feature/csv2shp.html).
Unexpectedly all the points given in the .csv file are being plotted at (0,0). I don't why is it happening.

Please help me resolve this issue.

I am using Maven and using the same code given in the Tutorial in the above mentioned URL.


texts: How do I load a .shp file using a button with ArcPy?

texts:
I have created a Python addin where the user should be able to select a shapefile and then load, but I've not been able to put in the current workspace.

import arcpy
import pythonaddins
import os
from arcpy import env as path
from arcpy import mapping

path = r"D:\MedellinPython"

class ButtonCarga(object):
"""Implementation for entregaFinal_addin.objButtonCarga (Button)"""
def __init__(self):
    self.enabled = True
    self.checked = False
def onClick(self):

    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
    addLayer = arcpy.mapping.Layer(path)
    arcpy.mapping.AddLayer(df, addLayer)

class ComboShp(object):
"""Implementation for entregaFinal_addin.objComboShp (ComboBox)"""
def __init__(self):
    self.value = "Seleccione archivo .shp"
    self.items = []
    self.editable = True
    self.enabled = True
    self.dropdownWidth = 'WWWWWWWWWWWW'
    self.width = 'WWWWWWWWWWWW'
def onSelChange(self, selection):
    pass
def onEditChange(self, text):
    pass
def onFocus(self, focused):
    objComboShp.enabled = True
    self.items = []
    for file in os.listdir(path):
        if file.endswith(".shp"):
            self.items.append(file)
def onEnter(self):
    pass
def refresh(self):
    pass


texts:
I'm trying to load raster file JPG format into QGIS on Linux. And then I got this message:
Invalid Data Source : /home/name/map/13-36.jpg is not a valid or recognized data source.

If I use "Add Raster Layer", then it say:
Invalid Layer: GDAL provider Cannot open GDAL dataset /home/name/map/113.jpg: libjpeg: Wrong JPEG library version: library is 62, caller expects 80 Raster layer Provider is not valid (provider: gdal, URI: /home/name/map/113.jpg

Same problem on QGIS 3.18. But on QGIS 3.10.4, It's complete normal.
Can anybody help me?
System info:
i7 4720 HQ
Intel Corporation 4th Gen Core Processor Integrated Graphics Controller
Linux mint 20.1
Cinnamon 4.8.6


texts: How to automatically reload raster layers if source is changed in QGIS?

texts:
How to  automatically reloaded raster layer if source is changed? (The path is the same but inside source file are changes)

I use:

 qgis.utils.iface.mapCanvas().refresh()

but layer not refresh.

#
#
#

The maps are repainting but there are still the same image (source file is changed on disc)

I use:

  layers = qgis.utils.iface.legendInterface().layers()
  for layer in layers:
      layer.triggerRepaint()


texts: What does `findCloseEdges` do in osm2po?

texts:
I have something like this:

int[] edges = graph.findCloseEdges(lat, lon, 1); // returns only one element
RoutingResultSegment rrs = graph.lookupSegment(edges[0]);
LatLon[] latlons = rrs.getLatLons()

For the red dot in the image (with lat/lon), array latlons contains the blue points. I would instead expect "close edge" to be a segment of the A3/E70 road.

Here is the image (I don't have enough reputation to post images apparently).

What is the third integer argument in findCloseGraph and what does this method return? I need to do simple "snap to road".


texts: Transform XY coordinates from within a custom flat geometry to latitude and longitude

texts:
I am looking for a Python based tool/API (PyProj?) that will help me transforming a set of XY points from within a custom flat map to latitude and longitude.

I have information about the latitude and longitude of (0, 0) point and also latitude and longitude of (Xmax, Ymax) point. The total area is not wider or longer than 400 meters.

It looks to me that PROJ/PyProj can do the job, however I am still lacking understanding how to initialize the thing.

Also, these XY points can be given as a matrix/table, so it would be nice if the tool can transform the whole data in vector implementation without the need to loop over all values.


texts:
In QGIS 3.8.2 I would like to split line into equal segments.
The issue is, that v.split splits the line into the required segments, except the last segment.

I would like to split the line into equal segments, but with the first and the last segment to be equal length. Any ideas?


texts: Mapinfo with Postgis - not installed in schema public

texts:
I'm aware of this

https://postgis.net/2017/11/07/tip-move-postgis-schema/

and I followed the steps and now I can see the tables in QGis; the problem now is MapInfo: I can open the table but I can't make it mappable with PostGIS (the 'postgis' option does not appear in the make table mappable option); is there some documentation about this?

Edit: There is a mapinfo.mapcatalog table in the database, the problem started when I changed the schema for postgis (from schema 'public' to another schema)


texts: Why doesn't Gdaladdo make the background transparent?

texts:
I noticed that the backgrounds of my geotiffs were transparent up until I started added overviews with gdaladdo.

Is there a way to combat this?`


texts: Which algorithm can I use to detect polygons from a set of line segments?

texts:
Say, I have a set of line segments. I need to check to see if they intersect to form any polygon or not. If they form any polygon, I need to detect and return them as a list.
Say, I represent a line segment as follows:
class Line
{
    public Point Start {get;set;}
    public Point End {get;set;}
}

and, say, a polygon is represented as follows:
class Polygon
{
    public List<Point> Points {get;set;}
    // ... ...
}

Which algorithm is generally used to solve this problem?


texts: loading a CSV with an embedded GeoJSON polygon geometry in QGIS

texts:
This open data portal offers download as CSV file. For larger datasets, that's the only way to download the whole dataset in one go. The CSV format is pretty creative: the polygon geometry is in GeoJSON in the "geometry" column.
If I drag and drop this file to my QGIS 2.18.1, it does not render as something geometrical. If however, I change the name of the "geometry" column to "WKT", QGIS does in fact render it. It only works with the WKT header, which is ironical. If I try to add the same info with the "Add CSV layer", it does not recognize the geometry column, as you can only feed that dialog window with WKT format, not GeoJSON.
Is this somehow normal behavior or is it a QGIS bug? Is it a valid option to offer this kind of format to the world on an open data website, or should I raise an issue to change the geometry data to WKT format?
Direct link to example data. Data structure:

"id";"objectid";"geometry";"gisid";"naam";"shape";"shape_length";"shape_area"
"1";"800";"{""type"":""Polygon"",""coordinates"":[[[4.4322857553331,51.236057528816],...[4.4322857553331,51.236057528816]]]}";"LEZ_01";"Lage-emissiezone (LEZ)";"";"1015.8237876604";"29420.558095414"


texts: Drawing polygons and quantifying them in buffer zone?

texts:
I am relatively new to GIS.

I need to define small areas (manually draw them out) within a buffer zone around a point. After they are drawn out I need to quantify how much of the buffer zone (1 mile from point) these drawn areas take up in total.

I have a ortho high-res map layer as well as a shapefile from the census. I also have several points geocoded for this project.

Can anyone help with the process of:

using the edit feature to draw out areas around the geocoded point. I have been able to create polygons but not sure how to create a feature that merges all drawn polygons that I am drawing as one feature.
seeing how much area these manually drawn areas take up in a buffer zone around the geocoded point.


texts: .mbtiles file with PBF to .mbtiles with PNG or JPEG

texts: I need this conversion to embed the mbtiles file in a mobile app. I know there are server based rendering solutions, but I just need to make this conversion once for a single file with OpenStreetMap .pbf vectors. What is the easiest solution?


texts: Layer Swipe OpenLayers 5 - How to change the left side?

texts:
I used the Layer Swipe for an app with OL5, like this example and it work perfectly for compare the baselayer with a new layer on the right side. But how to set a layer for the left side to ? I could change baselayer, but it appear behind the layer on right side, and I don't want.

var anneeCouche = "2016_20cm";

var _rasterSource = new ol.source.TileWMS({
  url: 'http://serveradress/geoserver/raster/wms?',
  params: {LAYERS: 'raster:ortho'+anneeCouche}
});

function changeDateRight(annee) {
  if (annee === '2013') {
    var anneeCouche = "2013_20cm";
  } else {
    var anneeCouche = "2016_20cm";
  };
  _rasterSource.updateParams({LAYERS: 'raster:ortho'+anneeCouche});
};

var _rightLayer = new ol.layer.Tile({
  title: 'Ortho',
  source: _rasterSource
});

_rightLayer.on('precompose', function(event) {
    var ctx = event.context;
    var width = ctx.canvas.width * (swipe.value / 100);
    ctx.save();
    ctx.beginPath();
    ctx.rect(width, 0, ctx.canvas.width - width, ctx.canvas.height);
    ctx.clip();
});

_rightLayer.on('postcompose', function(event) {
    var ctx = event.context;
    ctx.restore();
});


texts:
I have a set of files that have been mosaiced into a global VRT in a particular projection and resolution. However, I have a need to report statistics on a (say) one degree grid, where the reported statistics will be e.g. count, mean, mode or some percentile of the input pixels that fall within the grid cell. To make things more fun, there are NODATA instances too.

My solution is to create the mosaiced VRT in WGS84 Lat/Long, enforcing an output pixel size so that 1.0 % output_pixel_size == 0 (in other words, I can get an integer number of pixels into a 1 degree cell). Then, you run gdaladdo on the VRT and extract what you want. This should work for averages, but percentiles and some other things are a bit more complicated.

I could do it in Python, but I was wondering if there's an off-the-shelf solution I could use?


texts: Finding out if coordinate is within shapefile (.shp) using pyshp?

texts:
I'm trying to map coordinate data to a municipality. Reverse geocoders aren't as accurate as I need them to be (I was using geocoder pkg before), so I found a .shp that maps and names all the municipalities I need to refer to.

I have 9000 records in a CSV that need to be mapped, and was wondering If anyone knew how to check if the individual coordinates are in these shapes.

I'm not familiar with .shp files, and don't know if this would be possible through a package like pyshp

With help, I got it. Here is the code in case anyone has a similar issue in the future, it's sloppy and I may update when I'm done cleaning it up. Loading up the library of shapes is pretty slow, just to warn. And checking each point against each shape has high complexity, but for my sample of 9000 values against ~45 shapes it ran pretty quickly:

import shapefile as shp
from shapely.geometry import Point
import csv
from shapely.geometry.polygon import Polygon

fileName = 'TestLatLong'                #Open CSV file with points to check
sf = shp.Reader("NewMuni/NewMuni111")   #Open Shapefile with shapes to check points against
sfRec = sf.records() #Read records in shapefile

n = 0
m = 1
coor = ''
coorDict = {}
matplotDict = []
muniFinal = {}

for shape in sf.shapeRecords(): #Iterate through shapes in shapefile
    x = [i[0] for i in shape.shape.points[:]] #Initially for use in matplotlib to check shapefile
    y = [i[1] for i in shape.shape.points[:]] #Initially for use in matplotlib to check shapefile
    for i in x:
        matplotDict.append((x[x.index(i)],y[x.index(i)])) #Convert coordinates to be read by Shapely pkg

    munishp = Polygon(matplotDict)
    muniFinal[sfRec[n][1]] = munishp #Store shape in dictionary with key of municipality

    matplotDict = [] #refresh coordinate store for next shape
    n += 1

n = 0

with open(fileName + '.csv') as csvfile:
    reader = csv.DictReader(csvfile)
    for row in reader:

        coor = (row['latitude'],row['longitude'])
        rlat = float(row['latitude'])
        rlong = float(row['longitude'])
        if coor == ' ,   ' or coor == ', ':
            coorDict[row['PrimaryKey']] = 'No Data' #PrimaryKey is my primary key that I will use to write the data back into the .csv
        else:
            if float(row['longitude']) > 0:
                coorDict[row['PrimaryKey']] = (rlat,rlong)
            else:
                coorDict[row['PrimaryKey']] = (rlong,rlat)
        m += 1

#proof of concept- save the results however you'd like
for j in coorDict:
    for k in muniFinal:
        if muniFinal[k].contains(Point(coorDict[j])):
            print(j, 'in', k)


texts: Split vector layer in PostgreSQL/PostGIS

texts:
QGIS has this nice tool to Split a layer into multiple layers based on one attribute and create a new layer for each unique value within that attribute.

Is this also possible using PostgreSQL/PostGIS?
What would the SQL code look like?

Edit based on comments:
I do understand this goes against the idea of a RDBMS. Creating views is indeed a more handsome way.

I do need the multiple layers for one single reason. I am building a webmap where the legend should show the individual names of the polygons that are part of a set (the all have the same value for atttribute1, consider regions). For example: The original layer has the features 1,2,3,4, etc. with attribute1 is A, B, C or D, etc.

The new layers consists only of those features that have the same value in attribute1, but I do not want these features to be dissolved into one new feature/polygon. Then I put all those separate layers in a new map and style each layer with 'categorize on attribute2' but all with the same color so each feature within a layer all have the same color and each layer has a different color. Then the legend of that map looks like this:

Layer A
- 1
- 2
- 4
- 7
Layer B
- 3
- 5
- 6

I can do this within QGIS by using Spilt Vector Layer, and QGIS the creates six new files for each layer. Doing this with views in a database would be more convenient.


texts: Reading 13 bands (Sentinel 2A) .tif image in cv::Mat (opencv C++)

texts: I want to read 13 bands of a .tif image in separate cv::Mat objects for processing. I have tried reading it with opencv imread() function but it fails to load it. Also, I have used libtiff but I am unable to figure out how to get these 13 bands of images out of it in C++ for processing using opencv.


texts: ArcGIS Online Map App -- Table Widget -- Export Limit

texts:
Is there a way to change the Table Widget export limit from 2000 records?
I've found this:
https://support.esri.com/en/technical-article/000012383
But it doesn't do anything to change the limits of the Table Widget.


texts: Animate a series of polygons and points as a whole

texts: I have three KML files and each KML file contains a combination of points and polygons. I tried to animate this using timestamp but when I press play on Google Earth, it goes to individual points and polygons. I would like an animation that shows everything as a whole for each KML and not go through each point and polygon individually (meaning from a top view like the figure below). Basically I want to show how the polygons and points move across earth. Image below shows the three KMLs.


texts: Using SQLite instead of File Geodatabase?

texts:
I use ArcGIS and QGIS and I want to have common protocol and format of my files. To use a separate file for each program is not the best way, I don’t want to make a mess of data. I would ask if SQLite can play the role of gdb file. It is a format that can behave both programs the same (theoretically). I understand that there are differences with the gdb and SQLite but I need to manipulate in both software’s with only a file for local and online use.

Can I use SQLite instead of gdb?


texts: Change Table Field Editable to "yes" on a Hosted Feature Layer

texts: I have a hosted feature layer with one attribute field (text) having a 'no' within the settings editable property. I need to change this property to 'yes' within ArcGIS Online. Is this possible? Currently, this field has a lock symbol next to the field name.


texts: SPML (spatial panel models): Unbalanced panels

texts:
I am making spatial panel models from the splm library but I cannot make a fixed or random effects model because it comes out that it is an unbalanced panel.
Data
data and shp files
Reproducible example
library(dplyr)
library(Hmisc)
library(sf)
library(plm)
library(splm)
library(spdep)
data = read.csv("data.csv", colClasses=c(id="factor"))
map = st_read("nxparroquias.shp")
colnames(map)[1] = "id"

balanced = data %>%
  group_by(id) %>%
  #drop_na() %>%
  filter(all(c(2013:2018) %in% year)) %>%
  ungroup() %>%
  filter(id %in% unique(map$id))

map = map %>%
  filter(id %in% unique(balanced$id))
pt = spdep::poly2nb(map)
ptw = spdep::nb2listw(pt, zero.policy = TRUE)

describe(balanced$id)
#balanced$id
#n  missing distinct
#5832        0      972

describe(map$id)
#map$id
#n  missing distinct
#972        0      972

#checks if the data are balanced
plm::is.pbalanced(balanced$id,balanced$year)
#[1] TRUE

p_balanced = pdata.frame(balanced, index = c("id", "year"))  #same issues

Pooling
md1 = splm::spml(data = balanced,
                 y ~ x1+x2+x3,
                 index =  c("id", "year"),
                 zero.policy = TRUE,
                 model = "pooling",
                 listw = ptw)
summary(md1)

#ML panel with , spatial error correlation
#
#Call:
#  spreml(formula = formula, data = data, index = index, w = listw2mat(listw),
#         w2 = listw2mat(listw2), lag = lag, errors = errors, cl = cl,
#         zero.policy = TRUE)
#
#Residuals:
#  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#-14.860  -3.213   0.539   0.003   3.194  13.181
#
#Error variance parameters:
#  Estimate Std. Error t-value Pr(>|t|)
#rho 0.087841   0.019622  4.4767 7.58e-06 ***
#
#  Coefficients:
#  Estimate Std. Error t-value  Pr(>|t|)
#(Intercept)  5.1373431  0.0972298 52.8371 < 2.2e-16 ***
#  x1           0.0423977  0.0098276  4.3141 1.602e-05 ***
#  x2           0.1109780  0.0083800 13.2432 < 2.2e-16 ***
#  x3          -0.0444174  0.0096522 -4.6018 4.189e-06 ***
#  ---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

FE
md2 = splm::spml(data = balanced,
                 y ~ x1+x2+x3,
                 index =  c("id", "year"),
                 zero.policy = TRUE,
                 model = "within",
                 listw = ptw)

With reproducible example Fixed effects

Error in lag.listw(listw, u, zero.policy = zero.policy) :
object lengths differ

In my original data

Error in spfeml(formula = formula, data = data, index = index, listw = listw, :
Estimation method unavailable for unbalanced panels

RE
with random effects

Error in .C64("aplsb1", SIGNATURE = c(SS$signature, SS$signature, "double", :
NAs in argument 7 and 'NAOK = FALSE' (dotCall64)


texts: How do I create a SpatialLinesDataFrame from a dataframe?

texts:
I have created a data.frame in R . I want to get SpatialLinesDataFrame as a result so I understand that I have to convert the data.frame into lines, the lines into SpatialLines and the SpatialLines into SpatialLinesDataFrame. I did this

filedata.frame=data.frame(matrix(file),ncol=14, byrow=T))
file=Lines(filedata.frame)

I get the data.frame I want but not the lines. I am sure is basic but why I can't get the lines?


texts: Add buffer around polygon in meters using GeoDjango

texts:
I'm using GeoDjango with PostGIS. I have a table with a polygon field called "geofence".

I want to add a 5m buffer zone around the polygon.

obj.geofence.buffer(5)

This results in a huge buffer zone though. Experimentally it looks like the buffer zone is around 500km wide.

It appears that I must somehow do a transformation of my meter unit to the system used by polygon.buffer. How can I do that?


texts: How to import large ASCII into GIS (Grav/Mag) data?

texts:
I need to import two ASCII files (one for gravity data and one for magnetic data) into ArcMAP, but I can´t find any tutorial on how to do this, especially when the file is so large that I cant open it in either excel nor notepad? the file I can open in notepad contains X, Y, Longitude, Latitude, Elevation, FAA and BOUG.
I have tried to clean it and import, but I must be doing something wrong, what is the simplest way to import an ASCII file?


texts: How to make x and y tick labels to be printed as geographic coordinates in levelplot in R?


texts: How to modify contiguity matrix in ArcGIS?

texts:
Is there any way to modify the contiguity relationships between the polygons in ArcGIS. I am thinking about a topology database that lets the user modify the data on the map? For example, if there is no contiguity between an island and the mainland, you can change it on a topological level.

I have used an AZTool converter which generates an .AAT file which shows the relationships between the polygons in the data and lets the user modify those relationships. However, I do not think it's possible to reload the file to ArcGIS.

.AAT file showing the contiguity between the polygons


texts: Merge Raster Layers with Raster Calculator

texts:
I am looking to merge three raster layers together, we are only allowed to use the Raster Calculator tool in Spatial Analyst, ArcToolBox.

Using ArcMap 10.2


texts: How to change the thickness of this outline in QGIS?


texts: Example of a Dynamic Map Service that frequently changes?

texts:
A requirement of our web application is to automatically refresh a particular dynamic map service on a timer.

Requirement relates to a client who has a dynamic MS that consistently gets updates from GPS sensors.

Lets forget about Feature Services for now, and WMS or Image services that could change frequently.

Coding wise it was fairly straightforward to achieve this.

Our testers now want a Dynamic Map Service that changes frequently (ideally at least every 15s) so that they can see how this looks back in the front end.  Will also help us demonstrate the new capability to our client.

I have been scouring the Esri sample servers but cannot find one.

Is anyone aware of such a public service that we could use for some quick testing?

If the answer is no, there must be scripts out there to emulate GPS tracking movement, to perhaps keep changing the location of some points.  Then weI can setup my own service.


texts: Getting layers created in custom toolbox python script to add to mxd?

texts:
I'm trying to create a simple search tool so that my colleagues can easily query a database. My first version is like this.

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
SearchTerm = arcpy.GetParameterAsText(0)
Fuzzy = arcpy.GetParameterAsText(1)

if Fuzzy == False:
     NewLayer = arcpy.MakeFeatureLayer_management(r"W:\Source.lyr","Search Results (" + SearchTerm + ")","PROC_NAME LIKE '%" + SearchTerm + "%'")
else:
    NewLayer = arcpy.MakeFeatureLayer_management(r"W:\Source.lyr","Search Results (" + SearchTerm + ")","PROC_NAME = '" + SearchTerm + "'")

This works as expected if I just run the code through the python window (with the parameters hardcoded) but if I run it as a script in toolbox nothing happens. I suppose ArcMap doesnt automatically add the layer you create in a script, so I've added:

arcpy.mapping.AddLayer(df,NewLayer,"TOP")

This doesn't give add a layer to my mxd and gives me the error message: Runtime error <type 'exceptions.AssertionError'>:

Surely this should be pretty simple. Is there a better way of doing it?


texts: Problem in remove interaction after draw end in openlayers-3

texts:
I'm using OpenLayers-3.8.2. I have a draw Ineraction as follow:

function drawEnd(){
   map.removeInteraction(drawIneraction).
}
var drawInteraction = new ol.interaction.Draw({
   source: mysource,
   type: 'LineString'
});
drawInteraction.on('drawend', drawEnd);
map.addInteraction(drawInteraction);

I want to user draw just one feature. This Code work, but when user draw a feature, map zoom in 2 level.

How can I prevent zooming in?


texts:
To start off with:

Coding in C++

Using the QT Creator

Reading from WFS Version 1.1.0

Using Gdal/OGR compiled with Xerces

I'm trying to get all of the Capabilities from a WFS. I know how to do this via a Get Request, now I'm wondering how to do this via code (C++). So far I have worked out how to get the Feature Name and Feature Element ( Countries - String , CountriesFID - INT) etc... Now through code I want to be able to get the creator of the WFS, its version, Spatial references. And well basically everything from the getRequest(GetCapabilities) into object form for use in my program.

I'm currently using the Gdal "Hidden" layer called "WFSGetCapabilities"  which just returns the whole XML/GML document string form (not a list!) I would rather not have to QSplitString, I want my program to by dynamic and work out what's in any WFS and not just hard coded for one.


texts: Finding index of raster band with minimum or maximum value using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop and I have a 6-band raster.

How do I output a single-band raster whose values are the band number (1,2,3,4,5 or 6) of which band in the original raster contain the maximum (or minimum) value?

Accounting for ties would be a bonus.


texts: tool for point aggregation and mean calculation

texts:
I have the huge number of points spread across some map.

Is it possible to calculate the mean value of a parameter in the points of a small area and make one point for the same so that the total number of points in the map can be reduced?


texts:
I have a problem with aggregation of geometry either with ST_Collect, ST_Union, Array_Agg after doing intersection with transformed geometry.
Faulty PSQL code goes like this:

SELECT ST_AsGeoJSON(st_collect(st_transform(pos,94326))) as pos,array_agg(cnt) as total
from planet_osm_te_lvl_3
where ST_Intersects(pos,ST_Transform(ST_GeometryFromText('srid=94326;POLYGON ((14.0 44.0,  19.0 44.0, 19.0 46.0,  14.0 46.0,14.0 44.0))'),900913));

Column total is OK, but column pos returns null value. I figured out so far that the cause is the ST_Transform inside of ST_Intersect. The following code returns multipoint as expected:

SELECT ST_AsGeoJSON(st_collect(st_transform(pos,94326))), array_agg(cnt)
from planet_osm_te_lvl_3
where ST_Intersects(pos,ST_GeomFromText('POLYGON ((1737578.4691346379 5736562.273479164, 1737578.4691346379 5893061.453872256, 1932016.8447745405 5893061.453872256, 1932016.8447745405 5736562.273479164, 1737578.4691346379 5736562.273479164))',900913));

I have been busting my head over what could be the cause, but I can't figure it out. Any help would be greatly appreciated :)

OK, apparently there is possibly a buffer limit on the data that can be stored in a buffer/array. In my example there was around 500+ points that would return empty collection. Less than that would return multipoint.

https://stackoverflow.com/questions/9604832/how-to-convert-many-points-to-a-multipoint-in-postgis

Now I'm off to find a workaround :\


texts:
I was wondering if there is a way to draw a polyline but with an arrow head at the end point?

IPolyline line = new PolylineClass();
((IPointCollection)line).AddPoint((IPoint)sourceFeature.Shape);
((IPointCollection)line).AddPoint((IPoint)sinkFeature.Shape);


texts:
I have county shapefile and school district shapefile.

Once each school district file is overlay on county shapefile I would like to know how much percentage of area of school district file covered by each county shapefile.

Any idea how I can use GeoTools Java library or any other Java library to do that?


texts: How to change VB syntax to Python in Field Calculator?

texts:
I want to calculate a field, this is my working example in VB:

"C:\Users\Admin\Desktop\Data\" & !Work! & "\1000_\" & !Data! & ".pdf"

I'm slowly changing all my scripts to Python, how would I write this in Python?

got it:

'C:\\Users\\Admin\\Desktop\\Data\\'!Work!'\\1000_\\'!Data!'.pdf'

:)


texts: Python ogr2ogr language character encoding

texts:
In python script i'm exporting .shp files to .csv/.gpx/.kml

import ogr, ogr2ogr
...
ogr2ogr.main(["","-f", "CSV", csvfile, shpfile, "-lco", "GEOMETRY=AS_WKT", "-lco", "SEPARATOR=SEMICOLON"])
ogr2ogr.main(["","-f", "GPX", "-dsco", "GPX_USE_EXTENSIONS=YES", gpxfile, shpfile])
ogr2ogr.main(["","-f", "KML", kmlfile, shpfile])

The files are exported, but the language(Latvian) specific string field characters are incorrect. Shape files have string field "LABEL" with a value i.e. "marķieris".

In case of .kml file it is turned into "markieris", and in case of .csv and .gpx it is turned into "maríieris". All exported files have utf-8 (without BOM) encoding.

So far i had no luck finding if i have to pass some kind of additional parameter, or do something else to fix this language issue. Scrip is executed in Windows.


texts: Geoplot/GeoPandas missalignment in Windows (no Anaconda)

texts:
I am trying to familiarize myself with creating geographic plots in Python using geopandas and geoplot. So I am working my way through tutorials and documentation.

Currently, I am looking at the user guide for geoplot, but the created plot and the underlying map are misaligned. I am using Windows 10, Python 3.7 (no Anaconda), and am running inside Jupyter.

Here is the code

%matplotlib inline
import geopandas as gpd
import geoplot as gplt
import geoplot.crs as gcrs

continental_usa_cities = gpd.read_file(gplt.datasets.get_path('usa_cities'))
continental_usa_cities = continental_usa_cities.query('STATE not in ["AK", "HI", "PR"]')
contiguous_usa = gpd.read_file(gplt.datasets.get_path('contiguous_usa'))

ax = gplt.webmap(contiguous_usa, projection=gcrs.WebMercator())
gplt.pointplot(continental_usa_cities, ax=ax)

I am getting two future warnings and an image that is misaligned (see screenshot below).

What am I doing wrong, and what I can look at to debug/correct the alignment?


texts: Why are point-files cut at tiles-borders?

texts:
I'm placing a SVG in Tilemill via a point with lat and long.

So far that works great but when the tile the point is located in ends, the image is not rendered.

Is this a bug? Can I do something about it?


texts: Extract multiple indices values to the Feature collection in Google Earth Engine and export it as CSV File

texts:
var data= L8TOA1.filterDate('2019-02-01', '2019-03-28')
             .filter(ee.Filter.equals('WRS_PATH', 146))
             .filter(ee.Filter.eq('WRS_ROW', 40))
             .sort('CLOUD_COVER').first();
print(data);
//Get cloud cover information
var cloud= data.get('CLOUD_COVER');
print('cloud_cover is',cloud);
//Cip with boundary
var Study_region = data.clipToCollection(aoi);
//Indice
var NDBI = ee.Image (Study_region.normalizedDifference(['B6', 'B5']).rename('NDBI'));
var NDVI = ee.Image (Study_region.normalizedDifference(['B5', 'B4']).rename('NDVI'));
var SI= NDBI.addBands(NDVI)
print("Stacked_Image",SI);
var feature = ee.FeatureCollection([(Urban),(Vegetation)]);
print("Feature_points",feature);
//Extract values to the multiple points in different features
var img = SI.select(['NDBI','NDVI'])
 .reduceRegions({
   collection: feature,
   reducer: ee.Reducer.toList()});
print("Feature properties",img);
Export.table.toDrive({
   collection: img,
   description :"Indcies",
   folder:"GEE",
   fileNamePrefix: "Indcies_multifeature",
   fileFormat:"CSV",
});

I want to extract the values of NDVI and NDBI for multiple points in different featurecollection classes like Urban and Vegetation and export in different rows and columns. From the above code, the urban and vegetation are exported as 2 features(2 rows) and indices value of each point in the feature are exported together.


texts:
Currently pretty new to webmapping, I have created a PostGIS Database with OS Mastermap (OSMM) Data for Herefordshire (UK County). I have accsess and used these layers within QGIS so the database appears to be setup correctly.

I have successfully connected to the database within Geoservers web admin and been able to add other small vector layers with no problem. The OSMM layer currently don't show up as an option when I try to add them. Any ideas? File size perhaps?

I should mention this is currently setup on localhost.


texts: How to arcpy.Describe() a layer when script is run from outside of ArcMap?

texts:
The following code will work if I run it from the ArcMap Python dialog, or from a GP tool, but not from PyScripter IDE:

import arcpy

mxd = arcpy.mapping.MapDocument(r"c:\temp\<my_mxd_name>.mxd")
layerList = arcpy.mapping.ListLayers(mxd)

for lyr in layerList:
  print "Layer = " + lyr.name

  #Get all fields in the layer
  fList = arcpy.Describe(lyr).fields
  for f in fList:
    print "Field = "  + f.baseName

When I run the exact same code from PyScripter, it fails on arcpy.Describe(lyr).fields with the error IOError: "<TOC_layer_name>" does not exist, even though it does exist in the referenced MXD table of contents.  When PyScripter fails it opens up the ESRI file _base.py to the describe function.  I am assuming there must be something different about describing a layer in ArcMap vs calling it directly from Python, but cannot figure out what it could be.

How do I properly use the arcpy.Describe() function on a layer when I am not running the script directly from ArcMap?

I have figured out that if I use arcpy.Describe(lyr.dataSource) it functions properly, but only returns the fields from the source data layer.  However, since I'm trying to read fields from a joined table (which does work in the original code above when running from AcrMap), I really need to see the fields from all joined tables, not just the data source.


texts: QGIS 2.18 Make population point proportional

texts: I am mapping population abundance and distribution of a bird species in a one area. I am using several years worth of data to show how population and distribution vary by year. The years are marked differently by color. How do I make the year's population radius proportional to the year with the highest population? I need to do this so that the radii of the points are normalized.


texts: Adding ArcMap standard toolbar to the form in C# .NET WinForms project

texts:
I have Windows Forms project with ArcObjects framework.

I have already added axMapControl1 to the form.

The question is How can I add a standard toolbar to the form and link it with axMapControl1 ?

Like this:

Is there any standard toolbar control in the ArcObjects or should make it by myself?


texts: Writing MakeFeatureLayer_management Where Clause With Multiple Mixed Values using ArcPy

texts:
I am trying to select certain points in a shape file using three different fields. Two fields don't need to vary, but the third one (the analog_years variable below) does. I want to be able to select different values and numbers of years for this variable.
analog_years = [1950, 1993, 2010]

arcpy.MakeFeatureLayer_management(hurr_p_data_01, "lfalls_01",'"BASIN" = \'NA\'' " AND " '"USA_RECORD" = \'L\'')

How do I also include the analog_years variable in the where clause or is there a way to write a where clause variable for all three and just reference that instead?


texts: Calculating and mapping gains and losses of mangrove extent over time (vector data) using QGIS

texts:
I have a shapefile containing polygons (area of mangrove) per year for 35 years (=35 shapefiles).
How do I map changes of mangrove extent over time?

I want to map gains and losses (colour coded) between 1984 and 2019.

I would like to illustrate changes over time like this map for example.

For overall change, I have been trying to union the mangrove_2019 with mangrove_1984, creating a field 'difference' and calculate the difference area. Though I can't figure out how to union layers based on a corresponding field/column. And then I want to somehow visualise the gains and losses - perhaps by categorising the positive and negative observed changes and then symbolising these. Though this is not exactly working for me.
For yearly changes, my only idea is to rasterise the dataset first. But I've tried to do this and my output raster is blank.


texts:
I have a question about adding data to attribute tables. I uploaded an Excel file into ArcMap (File -> Add data -> Add XY data). I chose the table from Excel that had the orientation data's latitude, longitude, strike and dip. From there, I made the data permanent by exporting the data into a shapefile. This shapefile was then added as a new layer to my map. From there, using the attribute table in the shapefile, I was able to plot data as points on my map using the latitude and longitude. Then, I changed the symbol from a dot to one representing strike and dip, added rotation to the symbols using the strike, and lastly, labeled the symbols using the dip. This process worked well and I had no problems.

Now, my problem is whenever I edit the attribute table and want to add more data to that table, the new data does not show up on my map and I don't understand why.

Do I need to add that data to the existing Excel document and then go through the same process?


texts: How to determine objects area (size) measurement error using satellite images?

texts:
I am working at scientific paper about measuring object's size (it's area) using satellite images. I got some progress in it and need to compare my results with other scientific works. I have found some papers in russian scientific journals but can't find such papers in english.

Have you seen scientific papers about size (area) measurement error using satellite images?

Here is a picture of what happens when trying to measure some object's area
using satellite images:


texts: Earth Engine invalid band id when exporting

texts:
I am trying to export an image to Assets after adding a band representing time using the following code:

//Function to add a timestamp as a new band//
var createTimeBand = function(image) {
  return image.addBands(image.metadata('system:time_start').divide(3.154e10));
};

var img = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_062011_20130719');
var img = createTimeBand(img);

Export.image.toAsset({
  image: img,
  scale: 30
});

but I am getting this error after a few seconds when I try and run the task:

There seem to be no problems exporting without the timestamp band, but it is not clear to me why the timestamp causes this error.


texts:
Is it possible to load external wms services into ArcMap and publish them in ArcGIS Server?

I have tried this and get the error that the layer type is not supported.


texts:
I am performing a simple calculation on a polygon known to be approximately 6226 km^2 in area. It is stored in a Geography (WGS84 SRID) column.

The query is:

select st_astext(col), st_area(col) area from table

and returns:

"POLYGON((-180 58.282525588539,-178.916399160189 57.4759784390599,-178.191728834624 58.5761461944577,-180 58.282525588539))" | 5807028547.33813

The area returned (5807028547.33813) appears to be mm^2 and not km^2? The documentation http://postgis.net/docs/ST_Area.html states "by default area is determined on a spheroid with units in square meters"

Is this a documentation error, or is the above correct and I'm fundamentally misunderstanding the functionality?


texts: Which CRS would you use to georeference historical world maps?

texts:
Which CRS would you recommend to georeference historical world maps like the ones in the London Geographical Institute - The Peoples Atlas - 1920, e.g.

The map displays a lat/lon grid which I would like to use to pick control points. I managed to get some reasonable results simply using WGS84 EPSG:4326 but wonder if there isn't a better option.


texts: How do I calculate the latitude of an end point given a start point, bearing, and longitude of the end point?

texts:
I have a line that crosses the 180.  Given the bearing from Point A to Point B, how can I determine the latitude of the point on the line, given the longitude of the end point to be -179.9999?

For example:

Point A = 36N, 142E
Point B = 56N, 163W

If the initial bearing is 43° 31′ 51″, and my 'new' endpoint is at longitude 179.99E, how do I find the latitude?


texts: QGIS - Legend Columns - Uneven distribution of legend items

texts:
I'm working on a legend in QGIS 3.2 but cannot figure out how to get a set of legend items to distribute across 3 columns, but with an uneven number of rows.

For example I have 9 items that display in the legend, when I set the column count to 2, it creates a column with 5 items and a column with 4, however when I set the column count to 3 it give me 3 even columns with three items each.

Due to the layout of my map I want the items to display in 3 columns, but with 4 rows of items in the first 2 columns, and a lone item in the 3rd column.

Can anyone advise me on how to obtain this?

Current Outcome

Desired Outcome