texts: I want to scale, rotate, and move a polygon in a shapefile in QGIS 3.6 to align it with some other features. There are Rotate and Move tools in the Advanced Digitizing Tools but no Scale tool. After searching and reading, I found the v.transform tool but it seems to scale based on the CRS origin and not the center of the polygon. So, when I scale using the v.transform tool, it moves the polygon a long distance away from its original position, and I then have to shift it back near its original position again using the v.transform tool or the Translate tool. The scaling I need to do is sort of an iterative process, and it becomes very tedious if I have to scale and then shift multiple times. Is there an easy way to scale a polygon around its center using the v.transform tool or is there another tool that will accomplish this?


texts: Calculating distance between every entry in Table A and nearest record in Table B using PostGIS

texts:
I am new to PostGIS.

How do I calculate the distance and direction between every site (approx 30,000 record) and the nearest entity in a gazetteer (approx 300,000 records)?

Ideally a separate table would be created that contains the site_id, entity_id, distance between, direction between (entity & site). I have looked at an answer supplied by Tobias Herrmann in late 2014 - paraphrased below.

This code was written to do the calculation of distance, but is a work in progress.

SELECT site.gid AS gid_1, gazetteer.gid AS gid_2, ST_Distance(site.geom, gazetteer.geom) AS mindist
FROM table site, table gazetteer WHERE site.gid != gazetteer.gid AND ST_Distance(site.geom, gazetteer.geom) != 0
ORDER BY ST_Distance(site.geom, gazetteeer.geom)
LIMIT 1;


texts:
In QGIS do I have a EPSG:4326 shapefile that shows the coordinates "-7148300 -3688500" and I need them in the format -64º -31º for further processing in a geodjango app.

However using ogrinfo I get this information from the shapefile:

Feature Count: 22
Extent: (-64.279736, -31.509785) - (-64.092643, -31.332019)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]
id: Integer (10.0)

How can I convert the values?
I found useful both solutions:

1) Getting the coordinates of a given feature id, within Geodjango shell, so I have a reference and work from there using a custom conversion function.

2) Know the mathematical equivalence between the two pairs of coordinates.

A snippet or idea in geodjango would be nice, but any hint will be appreciated!

Thank you in advance.


texts: Making a single report per each feature in a feature class

texts:
I have been trying to figure out how to make an individual report for each record in a feature class.

The feature class is single point features and I do not want the report tool kicking out all of the information for each point feature that was collected in the field in a single report.

Rather, I would like each point to have its own report with its own attachments.

I have searched around on how to accomplish this but nothing has been helpful.

For example. When the report template runs it goes through each record, adds the information and its attachments then moves on to the next record. So what happens is that there is multiple records in one PDF file. I only want one record.


texts: Determining elevation of given coordinate (point) of rasterized .tiff-file in QGIS?

texts:
Currently, I'm developing a plug-in with Python for QGIS, which allows me to perform Line-Of-Sight Calculations. For the calculations, I'm using a python-script which calls modules like "r.viewshed" from GRASS to create the viewshed-maps.

R.viewshed requires some raster-file as input, and that's why I'm using  a method called "r.in.gdal" in my python-script to rasterize the elevation-map (tiff.file). The Line-Of-Sight Calculations are working properly, but there's just one problem:

r.viewshed requires one parameter called "observer_elevation". As for now, I've "hardcoded" this parameter since I don't know the elevation (height) of the given coordinate.

Does someone know, how to retrieve the elevation from a given coordinate so that I can process it in my python-script?


texts: Creating column in attribute table from summarised fields using QGIS?

texts:
I want to create a new (third) column in my attribute table. The fields of the new column are supposed to be the summarised value of the two other columns' fields. For example field 1 of the column I want to create is suppose to have the value 901 (663+238) and so on (se pic)

What function in the field calculator should I use and how do I write the expression?

It mostly seems to summarise the two columns and not the separate fields.. I'm not looking to use GROUP STATS or similar, but the field calculator only. I'm working with QGIS 2.18.15 on a Mac.


texts: Distance Analysis from single point to multiple points on vector layers

texts:
Is it possible to do the following in QGIS:

Create a text file that plots the distance from a single point (stored in one layer) and multiple points (stored in another layer).

I have tried to use the Vector Distance Matrix tool but without success so far. Maybe this is not the right tool?

Please note that I have no programming experience and am reliant on tools and plugins!


texts: Unable to load HYCOM ocean currents WCS layers in QGIS

texts:
I'm attempting to load HYCOM ocean currents model WCS layers into QGIS.  QGIS is able to access the server and the available layers at "https://wcs.hycom.org/thredds/wcs/GLBy0.08/latest?version=1.0.0".  QGIS is able to connect and retrieve the available layers, but reports the following error when attempting to add the layer:
"Invalid Layer: WCS provider Cannot get test dataset. Raster layer Provider is not valid (provider: wcs, URI: cache=PreferNetwork&crs=OGC:CRS84&dpiMode=7&format=GeoTIFF_Float&identifier=water_u_bottom&time=2020-08-17T12:00:00Z&url=https://wcs.hycom.org/thredds/wcs/GLBy0.08/latest?version%3D1.0.0"
I've attempted multiple HYCOM datasets with varying layer formats with no success.


texts: ArcMap Framework - Scroll bars in Layout View went missing. How to bring them back?

texts:
I have encountered a minor yet very annoying layout view change. All of a sudden my scroll bars in Layout View bars have disappeared.  I also no longer have the small layout view, data view, refresh and pause buttons which are normally a part of the the scroll bars. I now have to go to the main menu to switch my view from layout to data. However the scroll bars appear in data view.
I probably toggled something off accidentally but where? There must be a simple way to bring them back. Anyone?


texts: Error message Time Manager plugin with QGIS

texts:
I am using QGIS version 2.14.0 (Essen) and have downloaded the latest version of the Time Manager plugin in order to do some spatial-temporal analysis.

I have imported a csv file into QGIS and made a join wih a shapefile in order to import all the time data on the shapefile. However, when I use this shapefile in Time Manager, it recognises the data as the start date appears in the player time frame, but the video and animation do not load. In the logs I see the folowing error message:

Exception in function <function refreshTimeRestrictions at 0x0ECA18F0>:'ascii' codec can't encode character u'\xe9' in position 27: ordinal not in range(128)

I am using time data in the followng format, which should be compatible with Time Manager:

YYYY-MM-DD

Sample data:

start_date          end_date            start_date_auto     end_date_auto
2016-04-01 00:00    2017-11-01 00:00    2016-04-01 00:00    2017-11-01 00:00
2016-01-01 00:00    2018-01-01 00:00    2016-01-01 00:00    2018-01-01 00:00
2015-03-09 00:00    2016-12-31 00:00    2015-03-09 00:00    2016-12-31 00:00

Anyone know where I am going wrong?


texts: How to fill a field in a table using ModelBuilder ToolBox

texts: How to fill a field (set a value) in a table using ModelBuilder ToolBox.


texts: Solar Calculation and Projection in ArcGIS Pro on all buildings?

texts:
I am using this example to my project:

http://solutions.arcgis.com/local-government/help/solar-radiation-maps/get-started/explore-solar-radiation/

I have a layer of 3D Buildings as shown in the image:

When I applied the Calculation of the Solar Radiation on these Buildings Roofs, I found that the calculation is shown ONLY on the right 3d Building, whereas the middle and the left buildings are not calculated. Like this image:

Then, I suggested to myself to select these two buildings ONLY  and process the calculation of solar radiation to detect where is the problem. I did, but the result was not shown. Then I zoomed to the resulted layer, and found that these two 3D buildings are on Canada map rather being in Saudi.

I really surprising why this happened!

All of the 3D Buildings are in one layer as I mentioned earlier.

Is there any suggestion or solution why this problem occur in my work?


texts: Calculate Latitude based on a mouse position on Mercator map

texts:
I have a similar question to the one found here, Mercator projection: formula for latitude based on pixels? The answer was very helpful but my question is just different enough that it doesn't quite solve it.

Problem: Determine the latitude based on the mouse position on a map.

The map is 1000 pixels tall
Pixel 1 = Equator
Pixel 1000 = 80 degrees latitude

Is there an equation that determines at pixel X (for example 400) the mouse is at latitude Y (maybe 30 degrees)?


texts: Arcgis 10.2.2 Python Addin - Add Layer to Layout View?

texts:
I've developed a python addin toolbar that adds a desired layer file (chosen from a button nested in a menu in the toolbar) to the table of contents in a project.  The functionality only works, however, when the project is in the data frame view -- clicking a layer in the menu does nothing if the project is in the layout view.  I've googled around, but haven't found a solution to this issue.  Anyone have any ideas?

Using: ArcGIS 10.2.2, Python 2.7

Working sample code to add layer in Data Frame View:

import arcpy
import pythonaddins

class FD_Bau_Fischaufstiege(object):
    """Implementation for Themenmanager_addin.button (Button)"""
    def __init__(self):
        self.enabled = True
##        self.checked = False
        self.mxd = arcpy.mapping.MapDocument('current')
    def onClick(self):
##        layer = r'\\msds.wv.de\dfsroot\HV\GROUP\Gis\DATEN\wvdat\Querbauwerke\Fischaufstieg.lyr'
        layer = r'G:\Gis\DATEN\wvdat\Querbauwerke\Fischaufstieg.lyr'
        activeDataFrame = self.mxd.activeView
        df = arcpy.mapping.ListDataFrames(self.mxd, activeDataFrame)[0]
        if arcpy.Exists(layer):
            layerToAdd = arcpy.mapping.Layer(layer)
            arcpy.mapping.AddLayer(df, layerToAdd, 'TOP')
        else:
            warningButton = pythonaddins.MessageBox("Die Datei ist nicht verfugbar.\nBitte kontaktieren Sie die GIS Abteilung.", "Datei nicht verfugbar", 0)
        pass

I have tried the following, but it still won't add the layer when I'm in the Layout View:

import arcpy
import pythonaddins

class ButtonClass1(object):
    """Implementation for TEST_addin.button (Button)"""
    def __init__(self):
        self.enabled = True
##        self.checked = False
        self.mxd = arcpy.mapping.MapDocument('current')
    def onClick(self):
##        layer = r'\\msds.wv.de\dfsroot\HV\GROUP\Gis\DATEN\wvdat\Querbauwerke\Fischaufstieg.lyr'
        layer = r'G:\Gis\DATEN\wvdat\Querbauwerke\Fischaufstieg.lyr'
        activeDataFrame = self.mxd.activeView
        df = arcpy.mapping.ListDataFrames(self.mxd, activeDataFrame)[0]
        if activeDataFrame == 'PAGE_LAYOUT':
            if arcpy.Exists(layer):
                layerToAdd = arcpy.mapping.Layer(layer)
                arcpy.mapping.AddLayer(df, layerToAdd, 'TOP')
            else:
                warningButton = pythonaddins.MessageBox("Die Datei ist nicht verfugbar.\nBitte kontaktieren Sie die GIS Abteilung.", "Datei nicht verfugbar", 0)
        else:
            if arcpy.Exists(layer):
                layerToAdd = arcpy.mapping.Layer(layer)
                arcpy.mapping.AddLayer(df, layerToAdd, 'TOP')
            else:
                warningButton = pythonaddins.MessageBox("Die Datei ist nicht verfugbar.\nBitte kontaktieren Sie die GIS Abteilung.", "Datei nicht verfugbar", 0)
        pass


texts: Copying attachments from one geodatabase to another

texts:
I have a file geodatabase that is a compilation of multiple file geodatabases with varying schemas.

Some feature classes in the original databases had attachments, however the attachment tables no longer exist in the new, combined geodatabase.

There is an "Old-ObjectID" field for features in the new geodatabase that represents the ObjectID of the feature as it was in the original geodatabases.

I would like to know if there is some way to re-relate the old attachment tables to the new ObjectIDs in the new geodatabase?


texts: How to convert the SRTM DEM orthometric to ellipsoidal heights using ArcGIS or others GIS?

texts:
Recently I got the SRTM DEM from Cgiar,which based on WGS84/EGM96.

and I merged the tiles to one tiff file. I want to get the WGS84 ellipsoidal heights.

According the steps from the Arcgis Help 10.2 "Converting from orthometric to ellipsoidal heights",I do the whole steps,but I found the tif file was not changed,the height was the orthometric height based on EGM96 geoid. I Want one tif format DEM with ellipsoidal heights based on WGS84 ellipsoid.

I don't know where the problem is?

The steps I followed :

Create a file geodatabase

Steps:

Click the Catalog button  on the Standard toolbar.
This opens the Catalog window.
In the Location text box, type C:\arcgis\ArcTutor\Raster and press ENTER.
Right-click the Raster folder and click New > Folder.
Name the folder Exercises.
Right-click the Exercises folder and click New > File Geodatabase.
Rename the new file geodatabase ImageGDB.

Set the default geodatabase

Right-click the ImageGDB geodatabase in the Catalog window and click Make Default Geodatabase.

Create a new mosaic dataset

Steps:

Right-click ImageGDB in the Catalog window and click New > Mosaic Dataset.
Type HMA in the Mosaic Dataset Name text box.
Click the Coordinate System browse button .
Expand Projected Coordinate Systems >, choose wgs84 world mercator.prj, then click OK.
Click OK on the Create Mosaic Dataset tool dialog box.

Add rasters to the mosaic dataset

Steps:

Right-click the SRTM mosaic dataset in the Catalog window and click Add Rasters.
This opens the Add Rasters To Mosaic Dataset tool.In the Raster Type list, choose Raster Dataset.
Click the drop-down arrow and click Workspace.
Click the Input browse button.
Navigate to C:\arcgis\ArcTutor\Raster\Data\SRTM_tif and click Add.
Check Update Overviews.
Click OK to run the tool.

Add the Arithmetic function

Steps:

Open the Catalog Window  and navigate to the mosaic dataset.
Right-click the mosaic dataset and click Properties.Click the Functions tab.
Right-click the Mosaic Function and click Insert > Arithmetic Function.

Input Raster 1 is identified as the current contents of the mosaic dataset and will not be changed.
Click the Input Raster 2 browse button and navigate to the  \pedata\geoid folder in the ArcGIS install location.Choose the WGS84.img and click Add.

Click the Operation drop-down arrow and click Plus.
Click OK to close the Raster Functions Properties dialog box.
Click OK to close the Mosaic Dataset Properties dialog box.


texts: How to import/export features between QGIS and Electric Office

texts:
Background
I am working on a project that is looking at different electrical distribution networks across the UK. A utility company agreed to give me access to their data via their web map hosted using Electric Office Web. Unfortunately this was far too slow to be of any practical use so I was given the shapefiles that powered the webmap.
The problem
Unfortunately these shapefiles have been provided in such a way that QGIS pnly recognises them as a series of points. There is some metadata included in the attribute fields, but nothing I can use to perform any kind of analysis on to refine the mess of points.
I asked what GIS tool they used and the told me "EO" which I took to mean Electric Office. I have no knowledge of how EO functions or how it interprets data, or how it exports to shapefiles. Frustratingly, neither did the technician that had provided me the shapefiles.
The Data
Since I cannot disclose any information about the data, the project or utility company due to confidentiality, I have tried to set out as much information as I can without breaching confidentiality.
For a particular data set "Substation_locations". There are two (possibly three) accompanying shapefiles:

Substation_Locations (Pink points) is a point file and includes information such as the substation name, its geospatial location, ID, Asset Number, Asset ID, object ID but no data about the voltage of the network it is part of.
Substation_extents (Purple Polygons) includes the same information but provides a polygon for some of the substation boundaries. The ID, Name and Asset ID do match for the locations that have an accompanying polygon. Not much use as there is no additional data. There is also no data for substations that have no polygon.
Substation_annotations (Green points) is another point shapefile. Each substation_locations point is accompanied by three points in the annotation file. The attribute table contains fields; ID, Type, scale, text_size, Object ID. But the ID and Object ID do not match their respective substation ID's. Aside from their proximity to the substation_locations points there is no link between these two files and appears useless to QGIS.
connector_point_location (Yellow Points) another point file that seems to show locations where network connections are made between circuits. It has fields ID and Object ID, but again, these do not match any of the ID's in their respective substation_locations feature attributes. However, this file includes the information that is of value being the voltage of the connection. This data is not exclusive to substations and includes data for every connection between poles, towers, cabinets, meters etc. This makes it a more intensive process to perform Geoprocessing functions on.

I have included some images below to illustrate how this data appears within QGIS.
The Question

If you have knowledge of how EO works, is there a way to export data so that it can be displayed in QGIS as it would in EO? These files were supplied as shapefiles but QGIS can interpret quite a number of other file types. I can't seem to find any information regarding EO and its export capabilities.
Beside spatial joins between substation_locations and connector_point_location, is there a way to get the data in to the substation_locations attributes without creating many many duplicates or doing it manually?
I also have access to AutoCAD Map 3D, is this a better tool to use for this type of analysis?

substation_extent (purple) substation_location (pink):

substation_extent (purple) substation_location (pink) substation_annotations (green) connector_point_location (yellow):


texts: Create Vector Grid from a features Extent with PyQGIS

texts:
NOTE: THIS IS FOR QGIS 2.14.1 and uses VectorGrid instead of CreateGrid also though the question How to generate a grid programmatically in Python without GUI? is answered fantastically the PyQGIS language is outdated

We have a .txt file  "Extent.txt", in this example this file contains the coordinates of 4 points formatted like so;

  Point,Lat,Lon
  TL,-29.788254, 30.820720
  TR,-29.788254, 30.825372
  BL,-29.790710, 30.820720
  BR,-29.790710, 30.825372

We want to add this shapefile in our QGIS Map Document, and then create a vector grid that has the same extents of this .txt file, all using PyQGIS.
We add the .txt file with a PyQGIS script

import csv, ctypes, os, sys, datetime, string
import processing
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import *
from qgis.utils import iface

#Load CSV
gridSettings = QSettings()
gridoldValidation = gridSettings.value( "/Projections/defaultBehaviour")
gridSettings.setValue( "/Projections/defaultBehaviour", "useGlobal" )
gridInFlnm='Extent.txt'
gridInDrPth='PATH_TO_CSV_FOLDER'
gridInFlPth="file:///"+gridInDrPth+gridInFlnm
griduri = gridInFlPth+"?delimiter=%s&xField=%s&yField=%s" % (",","lon","lat")
gridbh = QgsVectorLayer(griduri, gridInFlnm, "delimitedtext")
gridbh.isValid()
gridbh.setCrs(QgsCoordinateReferenceSystem(32365,  QgsCoordinateReferenceSystem.EpsgCrsId))
QgsMapLayerRegistry.instance().addMapLayer(gridbh)
gridSettings.setValue( "/Projections/defaultBehaviour", gridoldValidation )

This produces our point vector layer from our .txt file!

So...How do we get our Vector Grid?

NOTE: THIS IS FOR QGIS 2.14.1 and uses VectorGrid instead of CreateGrid also though the question How to generate a grid programmatically in Python without GUI? is answered fantastically the PyQGIS language is outdated


texts: Reference a script/model from within the same Toolbox

texts:
in a script within a toolbox I'm trying to call a model in that same toolbox. Currently I have to hardwire the path.

arcpy.ImportToolbox(r"[PATH]", "LLPGToolbox")
arcpy.llpg_LLPGToolbox(inputFolder, output15, output21, output24, processedLLPG)

Is there a way of doing this without the need to hardwire paths each time (i.e. 'look in the toolbox I'm in for script X').


texts: Add dummy layer so that layer with "soft light" blending mode is visible by itself

texts:
I have a QGIS project with three raster layers:

SRTM GeoTiff (blending mode: SOFT LIGHT; color ramp: Viridis)
Google Terrain (blending mode normal)
Google Satellite (blending mode normal)

The problem is: if I display ONLY layer 1, nothing is visible, because it has nothing below to "blend with".

So my question is: is there a way to create a "solid color" layer, so that an upper layer with blending mode different than "normal" (and particularly with "soft light") can be visible somehow?

If I hide layer 2, I have the following result:

And if I display only layers 1 and 2 (layer 3 might be one or off), this is the result:


texts: I'm looking for a shapefile containing the locations of all the world's capitals as of (around) 2016?


texts: Easiest way to draw a lot of features over a base map (transparent tiles)

texts:
I'm working on a project, where I have to visualise a lot of lines (and  dots) on a map. We used to use openlayers, our own opensteetmaps tile server for "background" maps, and then generate a GeoJSON with 1-5k vectors and points in different colors, and it kinda worked fast enough in the browser (first 'render' took some time, but when zoomed in it wasnt too slow on a fast computer). It's hard to explain what exactly it showed, but imagine 'acceleration vectors' of a car driving down the road - so a point where the car was at the moment, and a line showing direction and acceleration force. Each car had a different vector/point color.

Then, things changed, and 1-5k became 1-5 million points+lines, and that number will probably still grow,  GeoJSONs are huge, and browsers just comes to a halt at 100% CPU usage.

So, how to draw them? I used to do something similar with points in php+gd and overlaying a google map, but it didn't show a vector if both (start, end) points were outside the visible field (but the vector crossed it). There could be better options for just that, but we don't know any.

The second option was making a tile server with transparent tiles and overlaying it over the osm maps, but except for following the tutorial for open street maps, we've never done that, and we don't know if it's a good/optimal idea.

I'd be also happy with any other ideas.

Our data is currently in a mysql database with DOUBLE lat/lon columns, but we have no issues with moving it to a postgis database.


texts: How to set Server path in GEOWEBCACHE_CACHE_DIR

texts:
I need to set server path in GEOWEBCACHE_CACHE_DIR inside Tomcat -> GeoWebCache -> web.xml.
The cache directory is on different server. I need to give IP/machine name and then the location.

How to set it? Please help.


texts: Map with custom tiles having uncommon EPSG

texts:
I've produced with mapnik my own tiles with EPSG:31256 on a shapefile.
I've structured the tiles in a OSM-like tree: zoom/x/y.png, with zoom from 0 to 2 and, x y starting from 1:

tiles
  ├── 0
  │   └── 1
  │       └── 1.png
  ├── 1
  │   ├── 2
  │   │   ├── 2.png
  │   │   └── 3.png
  │   └── 3
  │       ├── 2.png
  │       └── 3.png
  └── 2
      ├── 4
      │   ├── 4.png
      │   ├── 5.png
      │   ├── 6.png
      │   └── 7.png
      ├── 5
      │   ├── 4.png
      │   ├── 5.png
      │   ├── 6.png
      │   └── 7.png
      ├── 6
      │   ├── 4.png
      │   ├── 5.png
      │   ├── 6.png
      │   └── 7.png
      └── 7
          ├── 4.png
          ├── 5.png
          ├── 6.png
          └── 7.png

Tiles are 256x256 pixels and the tile at zoom 0 covers the entire shape, which is 7500x7500 meters (I used this dimension to set up scales). From zoom 1, tiles are splitted in 2 and so on. Map div element is 512x512 pixels.
To display tiles, I would like to use OpenLayers XYZ layer, with my own get_url function, as described here:

function create_map () {
    Proj4js.defs["EPSG:31256"] = "+proj=tmerc +lat_0=0 +lon_0=16.33333333333333 +k=1 +x_0=0 +y_0=-5000000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs "

    var proj = new OpenLayers.Projection("EPSG:31256");
    var bounds = OpenLayers.Bounds(5000,216250,12500,223750);
    var minZoom = 0;
    var maxZoom = 2;

    var map = new OpenLayers.Map({
        div: "map",
        allOverlays: true,
        projection: proj,
        maxExtent: bounds,
        allowSelection: true,
        scales: [7500, 3750, 1875],
        units: 'm',
        controls: [new OpenLayers.Control.PanZoomBar(),
                   new OpenLayers.Control.Navigation()]
    });

    var layer = new OpenLayers.Layer.XYZ(
        "my_tile_layer",
         '/tiles/',
        {
            maxZoomLevel: maxZoom,
            minZoomLevel: minZoom,
            getURL: get_my_url
        }
    );
    map.addLayer(layer);

    var center = bounds.getCenterLonLat();
    map.setCenter(center, minZoom);
};

function get_my_url (bounds) {
    var res = this.map.getResolution();
    // Note, I've added +1 because I do not have x and y coordinates equal to 0.
    var x = Math.round ((bounds.left - this.maxExtent.left) / (res * this.tileSize.w)+1);
    var y = Math.round ((this.maxExtent.top - bounds.top) / (res * this.tileSize.h)+1);
    var z = this.map.getZoom();

    var path = z + "/" + x + "/" + y + ".png";
    var url = this.url;
    if (url instanceof Array) {
        url = this.selectUrl(path, url);
    }
    return url + path;
}

My problem is that tile retrieval is wrong: at zoom 0, I correctly load 0/1/1.png, but it is positioned at the extreme right of the map: I had to pan it to make it visible. I do not understand why, because my center variable is correctly set to (8750, 220000).
If I zoom in I see nothing and I get not existent tile names such as 1/1/1.png.
Could you please help me?


texts:
Using QGIS 2.18.3, I have a project containing a single vector layer.  Both the project and the layer have the same CRS: NAD83/Alaska Albers (EPSG:3338).  I need to temporarily modify the project's CRS by changing the Longitude of Center from the default value of -154 to -163.

I've tried to edit the default longitude value in the Project Properties window, as highlighted in the screenshot below, but QGIS does not allow me to overwrite the default value.

Please note that I do not wish to reproject the data layer; rather I wish to modify the project display.


texts:
I'm wondering why after I updated to Python 2.7.6 on my OS X 10.9 Mac, QGIS Python console says that it is still running the old 2.7.5 version. This is not affecting my work but I'm curious about why this happens.

I though QGIS used the Python that comes with OS X but it seems it doesn't.

Do any of you know something about it?

I'm using KyngChaos QGIS 2.0.1 installers.


texts: Two versions of gdal_merge and two diffrent effects

texts:
Is there any possibility that two diffrent versions of gdal (1.7.3 and I believe 1.9.0) can give me two diffrent results with given:

two RGB files.Each twith 3 layers - input1.png and input2.png
command gdal_merge -o output_file.img -of HFA -separate input1.png input2.png

On version 1.9.0 i get 6 layered img file (3layers + 3 layers). On gdal 1.7.3 i get somehow 2 layers (not six) and i cannot recognize which layer is from which image.

These are really the same input files (i literally copied them to diffrent computer) and i used gdal from command line. What's funny is that it works good on windows, and on one ubuntu but on the other ubuntu it works not as predicted. It's the last idea i come with (diffrent gdal versions), but i find it very strange and unbelivable that same program works so much diffrently on diffrent versions.


texts: How can I login in GeoServer from a PHP application?

texts:
I'm looking for advice how to login on GeoServer from PHP.  I found much tutorial how to work with GeoServer through REST API and cURL: create workspaces, stores, etc. But still can't find how to login on GeoServer and still login to whole session. Is it possible?

In forums I found something but it does not work.

Any advice?

UPDATE

I change Acegi security to Spring Security and now code looks:

function loginGeoserver($username, $password) {

$geoserverURL = "http://localhost:8080/geoserver/j_acegi_security_check";

$post = http_build_query(array(
        "username" => $username,
        "password" => $password,
));

$context = stream_context_create(array("http"=>array(
    "method" => "POST",
    "header" => "Content-Type: application/x-www-form-urlencoded\r\n" .
            "Content-Length: ". strlen($post) . "\r\n",
    "content" => $post,
)));

$page = file_get_contents($geoserverURL, false, $context);

for($i = 0; $i < sizeof($http_response_header); $i++){

    $headerLine = $http_response_header[$i];

    $pos = strpos($headerLine, 'Set-Cookie');

    if ($pos === 0) {
            $str = explode("=",$headerLine);
            $value = explode(";",$str[1]);
            $cookieValue = $value[0];
            break;
    }

}

$cookieName = "JSESSIONID";
$cookieDomain = "http://localhost:8080";
$cookiePath = "/geoserver";
$cookieExpiration = 0;

setcookie($cookieName,$cookieValue,$cookieExpiration,$cookiePath);
echo $page;
return $cookieValue;
}

Now its return GeoServer index page but nothing more I still can't log in. How can I make Acegi authorization?


texts: What is getting nearest points between two different but overlaping datasets called?

texts:
I need some terminology help as I am still learning GIS.
I have two datasets.  One with the address of every fire and it's latitude/longitude of small census tract.  One with the location of every parcel of land (property) within that census tract.  I need to join or match the two datasets so I end up with a map of every property that has ever had a fire and a list of the properties without ever having a fire.  I thought it might be similar nearest neighbor or involve R-trees.  However,  every time I look for a way of doing this I run into methods of "matching point within the same dataset" not what I need.  I have to use two different datasets.
What exactly is this type of matching or joining between nearest lat/long of different datasets called in GIS terminology?
To help answer this question let me show two different dataset with slightly different situations:
(It took me a minute of searching to find these polygons)
First is the my actual GA data. This data is only of 3 months:

In it, every green point is a building fire, while every brown shape is a building.  About 95%-98% of the data overlays a parcel, so I think is a point in polygon overlay.


texts: Open Street Map Routing using Ships and Trains

texts:
I wanne know if is there a way to use the Open Street Map API to know the distances of ships and getting the distance of two points.

Maybe there is a profile for it i have searched for Marine Navigation and OSMR but i didn't found anything.

And Also can you install the Open Street Map Server on Windows.


texts: Converting GPX timestamps to different format for ArcGIS Desktop?

texts:
I'm trying to figure out an easy or efficient way to convert the timestamps contained in GPX files to a different timestamp. I have data for air quality levels saved in an excel spreadsheet and want to match up timestamps between my air quality data and my GPX file so that I can join tables in ArcMap and overlay air quality levels over the GPX track.

To give you an idea of what I mean, this is what the timestamp looks like in the GPX file:

  2017-11-07T12:38:55Z

And this is what the timestamp looks like as saved in the excel sheet for my air quality data:

  07/NOV/2017 12:38:55

If I could figure out a way to convert one timestamp to another without going through each row and changing them manually then I could use these timestamps to combine the tables from my GPX file and the excel sheet so that I can overlay my air quality data onto a map.

Is there a tool online, in ArcMap, or an easy way in excel to mass convert timestamps so that the two timestamps I posted above are in the same format?


texts:
GRASS has extensive and brilliant documentation, but I still have a bit of trouble using v.net.path. in the example for shortest path, it says: "

#Shortest path: ID as first number, then cat1 and cat2

echo "1 1 2" | v.net.path myroads_net_time alayer=3 nlayer=2 out=mypath"

I've had trouble understanding what ID is; my table has 'cat' as it's primary key, but there is no ID column.

Any ideas?

P.S. I'm using Grass 6.4 on Mac.


texts: Python QGIS plugin could not load on linux build in window

texts:
I have created a plugin in QGIS windows version using
"Plugin Builder". which on a button click opens another window, that is calling another class in the dialgoue python file, with another .ui file. That plugin is working fine on windows but when i installed that plugin on linux QGIS i got below error

Couldn't load plugin tplUserStatsReport due to an error when calling its initGui() method

TypeError: Expected bytes or unicode string, got newstr instead
Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 342, in startPlugin
    plugins[packageName].initGui()
  File "/home/hasan/.qgis2/python/plugins/tplUserStatsReport/tplUserStatsReport.py", line 203, in initGui
    self.setInitialStats();
  File "/home/hasan/.qgis2/python/plugins/tplUserStatsReport/tplUserStatsReport.py", line 330, in setInitialStats
    connect= psycopg2.connect("dbname='"+str(Configuration.dbname)+"' user='"+str(Configuration.username)+"' host='"+str(Configuration.host_ip)+"' port="+str(Configuration.port)+" password='"+str(Configuration.password)+"'")
  File "/usr/local/lib/python2.7/dist-packages/psycopg2/__init__.py", line 129, in connect
    dsn = _ext.make_dsn(dsn, **kwargs)
  File "/usr/local/lib/python2.7/dist-packages/psycopg2/extensions.py", line 155, in make_dsn
    parse_dsn(dsn)
TypeError: Expected bytes or unicode string, got newstr instead

Python version: 2.7.12 (default, Nov 19 2016, 06:48:10) [GCC 5.4.0 20160609]
QGIS version: 2.18.2 Las Palmas, 102ee87

Python Path:
/usr/share/qgis/python/plugins/processing
/usr/share/qgis/python
/home/hasan/.qgis2/python
/home/hasan/.qgis2/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python2.7
/usr/lib/python2.7/plat-x86_64-linux-gnu
/usr/lib/python2.7/lib-tk
/usr/lib/python2.7/lib-old
/usr/lib/python2.7/lib-dynload
/home/hasan/.local/lib/python2.7/site-packages
/usr/local/lib/python2.7/dist-packages
/usr/local/lib/python2.7/dist-packages/TileStache-1.50.1-py2.7.egg
/usr/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages/PILcompat
/usr/lib/python2.7/dist-packages/gtk-2.0
/home/hasan/.qgis2//python

How can i resolve this conflict that is occuring due to a plugin build on  windows and run on Linux?


texts: Assign values to points according to a pattern

texts:
I'm making a map of trees in QGIS (3.10). Each point represents a tree and the attribute table has a column that should contain the tree species. For now this column is empty, but I would like to fill this column. The trees are planted in a pattern, for example cashew tree - baobab tree - banana tree.

How can I assign the tree species to the table according to the pattern? So every first tree is a cashew, every second tree is a baobab and every 3rd tree is a banana. The fid of the table is in the right order, so I could build an expression with fid+0 = cashew, fid+1 = baobab and fid+2 is banana. Would that be possible?


texts: Adding consecutive unique ID starting from 1 for each distinct set of values

texts:
I need to add a unique ID (using $ID or @row_number) but the problem is that I need to start from 1 to "...n" considering identical values from other field, then, when there is a different value, start counting from 1 again, if the value is repeated, write 2, 3, 4, 5 until finding new distinct values, and so on.
I have an image more explanatory than me, the red column is what I am getting, the green column is how I would like to add the numbers. I guess using the expression "count_distinct" could be useful, but I don't know how to structure it with my values.


texts:
I have created a class library which does some geoprocessing. The addin calls a class that is an async process. I have ensured the thread is STA and the arcobjects are thread safe (ie. not passed from the addin). All arcobjects are created within the thread.

Because this is a class library I wrapped it in a winforms UI and also as an addin. Both sets of code are exactly the same and test were performed using the exact same data. Both the winforms and addin complete the code with the desired results and no memory leaks are evident. For the case of the addin there is no interaction with the map period at this point and there are no mapping or display elements in the winforms code either.

the only UI updates are the update of a progress dialog in both the addin and UI. The addin is using a dockable window (user control ui).

The issue I am seeing is when the library is called from the addin, code execution is 5x slower than the same code called through the winforms application.

Any ideas on where I might look to see why this is occurring?


texts: QGIS: when using composer to save as image the image has a shift of about 50m

texts:
I'm using QGIS 1.8.0; I have two layers turned on, one is the Google Earth hybrid layer from the OpenLayers plugin and the other a shapefile of my own making.

I set on the fly projection and everything looks spot on in the active window. When I got to composer everything still looks good. It's only once I have saved the image (exported) that there is an issue.

I think I read about this earlier but I can't find the answer to it anywhere. I can't remember if it has to do with on the fly projection or not... it seems ridiculous that the composer doesn't just "print the screen" of what's displayed there.

Any help would be great.

Shapefile is located here

On the fly projection is UTM 50S or EPSG:32750


texts: SQL 2012 code that spatially polyline fine (join file) to a point file (target) and provides Distance Field

texts:
I am currently working in ArcMap 10.2.2 but would like the SQL code to create the same output as that provided by the Spatial Join tool.

Target feature = point, Join Feature = polyline, Option CLOSEST and a Distance field generated and shown in new shapefile.

I have been struggling with the STDistance function but this is only for "a" distance between a single pair of geometries.


texts:
I would like to calculate the lenght(km) of each line in shapefile of roads in a study area. I tried to use "Calculate Geometry" but it doesn't work! I cannot select 'Units', and 'Property' for 'Length' is "Length - Disabled" in Calculate Geometry dialog.

How can I calculate length of roads?

Coordinate system is WGS84.


texts: all-in-one gis packages for windows 8.1?

texts:
Is there an all-in-one installer that allows having all open-source gis apps installed in windows 8.1? with all apps I mean for instance postgresql, postgis, mapserver, qgis, grass.

I have cygwin installed, but I've read that it's way better to install these apps in windows. In case I need to query my psql db located in my postgresql windows installation from cygwin (ie. shell terminal), could that be possible?


texts:
I have local tile also 20 level. I found out that OSM tile level 20 don't loading.

map code:

map = new ol.Map({
 target:'map',
 renderer:'canvas',
 view: new ol.View({
    projection: 'EPSG:3857',
    center:[2577400.0, 7007400.0],
    extent:[2567000.5, 7000400.5, 2588377.0, 7014361.0],
    minZoom: 13, maxZoom: 20, zoom:14 })]
  })
});

and my osmLayer code:

var osmLayer = new ol.layer.Tile({
       source: new ol.source.OSM({
       url: '../mapy/a_tiles/{z}/{x}/{y}.png'
 })
});

Is there any problem with OSM level 20?
In browser devTools I have:

Request URL:http://server/mapy/a_tiles/19/295869/170467.png
Request Method:GET
Status Code:200 OK

despite I just loading 20 level ...

Anybody have any solutions, any ideas?

20-level zoom view:

19-level zoom view:


texts: Error when EXPORTING Merged Raster with QGIS 3.14

texts:
I have been merging OS Open Data Terrain 50 asc files to create a single Geotiff file. I use Raster> Miscellaneous > Merge...
This works flawlessly if in the Merged slot [Save to temporary file] I choose Save to File... (Geotiff)
However, if after the merged layer has been generated, I EXPORT the layer to a Geotiff it creates a file with numerous grid square artifacts. This didn't use to happen. Has something changed with the export option?


texts: Update field based on spatial query (QGIS)

texts:
I was an intermediate user of MapInfo and have recently converted to QGIS. I would often use MapInfo to add a column in boundary file A saying which zone in boundary file B it is contained in. However, I cannot work out how to do this in QGIS.

Sorry, this is the error message I received when I tried to run the function.


texts:
I wanted to add ILayer to IMap, i think it has been added but problem is that it does not display it on the website.

I am running ArcGIS Server 10, using visual web developer 2010 express edition for web mapping application, vb.net (web adf). Here is the code:

 Public Sub ServerAction(ByVal args As ToolEventArgs) _
    Implements ESRI.ArcGIS.ADF.Web.UI.WebControls.Tools.IMapServerToolAction.ServerAction

    If TypeOf args.Control Is ESRI.ArcGIS.ADF.Web.UI.WebControls.Map Then

        Dim map1 As ESRI.ArcGIS.ADF.Web.UI.WebControls.Map
        map1 = CType(args.Control, ESRI.ArcGIS.ADF.Web.UI.WebControls.Map)

        Try
            Dim mapFunc As ESRI.ArcGIS.ADF.Web.DataSources.IMapFunctionality = _
                CType(map1.GetFunctionality(0), ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.MapFunctionality)
            Dim mrl As ESRI.ArcGIS.ADF.Web.DataSources.ArcGISServer.MapResourceLocal = mapFunc.Resource
            Dim sc1 As ESRI.ArcGIS.Server.IServerContext = mrl.ServerContextInfo.ServerContext
            Dim mapS As ESRI.ArcGIS.Carto.IMapServer = mrl.MapServer
            Dim mapSObj As IMapServerObjects = mapS
            Dim pMap As ESRI.ArcGIS.Carto.IMap = mapSObj.Map(mapS.DefaultMapName)

            MsgBox(pMap.Name & ": " & pMap.LayerCount)

            Dim f As ILayerFile = New LayerFile()
            f.Open("C:\data\a.lyr")
            Dim pLayer As ILayer = f.Layer
            pMap.AddLayer(pLayer)

            MsgBox(pMap.Name & ": " & pMap.LayerCount)

            mrl.RefreshServerObjects()
        Finally
            map1.Refresh()
        End Try

    End If
End Sub

Please can you explain or show why the added layer is not being displayed?


texts: pgrouting 2.0 on qgis cannot create route

texts:
I am trying to create a network using pgrouting and PostGIS.

I was initially using OS Open Roads data but even when I follow this tutorial (https://anitagraser.com/2013/07/06/pgrouting-2-0-for-windows-quick-guide/) and the same data set all i get is the following;

The issue is that when I now attempt to run the route nothing happens, I am using the suggested SQL;

SELECT seq, id1 AS node, id2 AS edge, cost, geom
  FROM pgr_dijkstra(
    'SELECT id, source, target, st_length(geom) as cost FROM network.publictransport',
    1, 3000, false, false
  ) as di
  JOIN network.publictransport pt
  ON di.id2 = pt.id ;

But, when I execute nothing happens, I do not get an error nor is it loading it just returns 0 columns instead of them being populated.


texts:
How can I retrieve database objects from postgis database which are contained inside Leaflet circle object? I assume circle is only a point with some radius and not a polygon. I wonder can it be somehow converted to the polygon so I can read objects inside it from postgis database.

Thanks!


texts: ArcGIS equivalence of QGIS autofill field using *Default value* by expression and spatial relations

texts:
I am trying to find an ArcGIS equivalent for the QGIS operation that allows to fill automatically a field during the digitising using the option Default value by expression in Proprieties > Attribute Form with one expression based on spatial relation between two layers.

In QGIS I am using the function aggregate in an expression that permits to autofill the field of a line when the geometry touches a point layer that contains the information.
This is the QGIS expression used in the image example:
array_to_string(
aggregate(
'point',
'array_agg',
"try 02",
intersects($geometry, geometry(@parent)))
)


texts: Finding which polygon feature is most in using PyQGIS?

texts:
I have a vector layer of land parcels and another layer of polygons separating them into territories. I am trying to identify which territory each parcel is in within a python plugin. My plan is to iterate over the parcels, and for each parcel iterate through the territories and test which it's inside. I'm not sure whether to use .contains(), .within(), or .intersects()

My concern is what happens if a parcel is half inside one territory, and half inside another. I'd like the plugin to figure out which has the biggest portion of it and return that one. How do I test for that?

I've thought about using a spatial join to accomplish this, but for some reason I get errors with my dataset. And I want more control in this exact situation than arbitrarily choosing the minimum or maximum value for object ID)


texts: Extract By Mask moving pixels in ArcGIS Desktop?

texts:
I am creating some statistics based on certain values. The issue I am having is that the extracted values are not the same as the values in the original raster. This creates incorrect output, and I need to ensure that, when joined, both fields are the exact same. After looking on ArcGIS Desktop, i found that the reason is that the pixels are moved by a tiny bit on the extract vs the original.

Hopefully you can see that the pixels are not aligned properly, and I need them to be properly aligned. I have tried the following:

Set extent and Co-ordinate system to the extraction layer (in Extract by Mask)
Snap the raster to the extraction extent(in Extract by Mask)

I am not sure what else I can do. I understand that this is a glitch as documented http://support.esri.com/en/knowledgebase/techarticles/detail/34448

From the article, setting the extent should fix it, but it does not do anything.


texts:
I have a marine AIS feed that is being decoded and written to a Google Earth KMZ file using aisDispatcher software. I am trying to get QGIS 3 to update the layer automatically, but am constantly failing. I have tried the rendering/fresh layer option, but this seems to do nothing.
I've tried a lot of python examples, but nothing seems to update or refresh. The only thing that works to refresh is F5. To that end, I found the below to refresh in Python console:
from qgis.utils import iface
iface.mapCanvas().refreshAllLayers()

However, it is slightly outside my comfort zone to edit this to put an interval timer in to put it into a loop.
Similar thread here:
Automatically refreshing mapCanvas


texts: How to convert the BLOB value for the features geometry of type st_geometry kept in ArcGIS databse on Oracle

texts:
I am using ArcGIS Oracle database and I am trying to access feature geometry points as string. But the problem is that the points are kept as blobs and I can't convert to meaningful string.

select  sde.ST_AsText(a.SHAPE)
from AKBIS_GEO.A236 a

About query returns response like below:

    ORA-28575: unable to open RPC connection to external procedure agent
ORA-06512: at "SDE.ST_GEOMETRY_SHAPELIB_PKG", line 67
ORA-06512: at "SDE.ST_GEOMETRY_OPERATORS", line 67
28575. 00000 -  "unable to open RPC connection to external procedure agent"
*Cause:    Initialization of a network connection to the extproc agent did
           not succeed. This problem can be caused by network problems,
           incorrect listener configuration, or incorrect transfer code.
*Action:   Check listener configuration in LISTENER.ORA and TNSNAMES.ORA, or
           check Oracle Names Server.

When we try this query:

select  a.SHAPE.points
from AKBIS_GEO.A236 a

The text value corresponding to the blobs related to points shows as follows:
(which is not useful for me)


texts: Distance between point and line layers & spatial index issue

texts:
I am trying to compute on QGIS the distance from some clusters (ie points) to their closest road (ie lines);  no need for all combinations, just the distance from each cluster to the closest road. To do so, I saw on a forum (on this post) that a good method was to transform lines into points, which I did, and then to compute the distance between both point layers. Nonetheless, when trying to do the last step, I get this error :

"/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/tools/vector.py", line 211, in spatialindex

idx = QgsSpatialIndex(layer.getFeatures(request))

Exception: unknown"

I am using a MAC. I tried on both QGIS 2.18.2 and QGIS 3.0.3, and got the same error. I tried, as advised on a page, to export geometric attributes using MMQGIS, and then re-importing them on QGIS, but same error. I also heard that it could be useful to check the validity of points ; nonetheless as I am using all the major roads network in Africa, I get millions of points, and checking them would take hours and hours ... I tried on QGIS3 to compute the distance between points and line centers (which does not give me the answer I want, as it is not the closest distance, but I did so to check if it was working this way), and it worked well, so I guess the issue lays in the transformation of my line layer into points

Has someone faced this issue of spatial index error, and knows what is means exactly? Is there a way to correct the points, to get rid of those not correct if this is the issue ?


texts: QGIS: Rotate SVG symbol in layout according with degree number in another layer

texts:
QGIS 2.18.3:

I have a shapefile layer in QGIS which describes points. One class of points are ship lane beacon lights (svg image), and these should be turned towards the ship lane. The angle of the ship lanes are found in another layer (the shapefile line layer for the ship lanes) and the angle is in degrees with decimals (for instance 248.8000 ie 248.8 degrees).

It is possible to edit the layout rule for the SVG symbol, and there I can edit the rotation calculation. But how do I tell QGIS to use a field "ANGLE" in another layer "Shiplane" for the rotation angle?


texts: How to extract drainage network through gdaldem

texts:
From gdaldem I can view slope, aspect, and hillshade maps.

"gdaldem slope manj.tif ouoi.tif -p -s 60"
"gdaldem aspect manj.tif gg.tif"
"gdaldem hillshade manj.tif hh.tif"

Is there any command to get flow direction map from gdaldem?

Can you give a working example of flow direction and drainage network extraction?


texts: Calculate Color of Object So No Adjacent Objects Repeat

texts:
I am attempting to assign one of 16 colors to map objects (e.g. counties, states, ZIP Codes, etc.) so that no two adjacent objects share the same color.

I'm seeing tons of articles online about "edge coloring" and the "four color theorem," but I cannot seem to find a way to apply these theorems as an algorithm in PostgreSQL/PostGIS (or any practical application, for that matter).

I'm 100% sure this is a solved problem, but since Google isn't revealing an answer, I'm guessing I'm too ignorant on this topic to ask the right question.

Can someone point me in the right direction?


texts:
I'm referring to this post:

Why does QGIS 3.2 "native:extractvertices" algorithm not work properly in standalone script?

as part on creating startup file py3-env.bat for opening PyCharm with PyQGIS works:
@ECHO OFF

set OSGEO4W_ROOT=C:\OSGeo4W64

@echo off
call "%OSGEO4W_ROOT%\bin\o4w_env.bat"
call "%OSGEO4W_ROOT%\bin\qt5_env.bat"
call "%OSGEO4W_ROOT%\bin\py3_env.bat"

@echo off
path %OSGEO4W_ROOT%\apps\qgis\bin;%PATH%
set GDAL_FILENAME_IS_UTF8=YES

set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins

SET PYCHARM="C:\Program Files\JetBrains\PyCharm Community Edition with Anaconda plugin 2020.1\bin\pycharm64.exe"

set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python37
set PYTHONPATH=%OSGEO4W_ROOT%\apps\Python37\lib\site-packages;%PYTHONPATH%

set QT_QPA_PLATFORM_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\Qt5\plugins\platforms
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT%\apps\qgis

cd /d %~dp0
::python3 scratch.py
::pause
start "PyCharm aware of QGIS" /B %PYCHARM% %*

I have:
QGIS 3.12.1 installed in: C:\OSGeo4W64
PyCharm 2020.1 installed in: C:\Program Files\JetBrains\PyCharm Community Edition with Anaconda plugin 2020.1\
After loading PyCharm the interpreter looks like this:

Now I just would like to run a simple test script like:
#native:extractvertices
import sys

from qgis.core import (
     QgsApplication,
     QgsProcessingFeedback,
     QgsVectorLayer
)
from qgis.analysis import QgsNativeAlgorithms

QgsApplication.setPrefixPath(r'C:\OSGeo4W64\apps\qgis', True)
qgs = QgsApplication([], False)
qgs.initQgis()

sys.path.append(r'C:\OSGeo4W64\apps\qgis\python\plugins')

import processing
from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

layer = QgsVectorLayer(r"D:\test_pygis\lines.shp", 'my layer', 'ogr')
output = r"D:\test_pygis\verticles.shp"
params = {
    'INPUT': layer,
    'OUTPUT': output,
}
feedback = QgsProcessingFeedback()
res = processing.run("native:extractvertices", params, feedback=feedback)
print(res)

First the processing module is not found:

If I try to run the code the second well know error appears in the Python console:
qt.qpa.plugin: Could not find the Qt platform plugin "windows" in "C:\OSGEO4~1\apps\Qt5\plugins\platforms"
This application failed to start because no Qt platform plugin could be initialized. Reinstalling the application may fix this problem.


texts: Getting Infographics within shapefile boundary using ArcGIS Pro

texts:
I just got ArcGIS Pro and have been googling around trying to figure this out.

How do i get the Infographics that are within a shapefile?

I would like to know the population that lives within the shapefile


texts:
I am working with ArcMap 10.2. In the following code, I have a function, build_dictionary, that builds a dictionary that will be used to populate a field within another layer. With the second function, update_points, I am trying to use an UpdateCursor to iterate through each feature and if that feautre's OBJECTID is a key within the dictionary, populate a field called "RoutingStreetExtKey" with that key's value.

I have verified that the build_dictionary function works properly with an "AddMessage" (and therefore omitted the code). I return the dictionary in order to use it with the next function.

The problem occurs with the update_points function. As the script runs, I can see it iterate through the point layer, however it only updates the last feature. All other features have nothing written to the "RoutingStreetExtKey" field except the very last point feature. The last feauture is populated with the correct value.

I am concerned that I may not be passing the dictionary correctly from build_dictionary to update_points. Furthermore, is there something else I need to add to the body of UpdateCursor to ensure all features are being updated?

import arcpy
point_layer = arcpy.GetParameterAsText(0)
street_layer = arcpy.GetParameterAsText(1)

def build_dictionary():
    #Code...
    arcpy.AddMessage(ext_key_dict)
    return ext_key_dict

def update_points(passed_dict):
    with arcpy.da.UpdateCursor(point_layer, ["OBJECTID", "RoutingStreetExtKey"]) as cursor_update:
        for new_row in cursor_update:
            obj_id = new_row[0]
            if obj_id in passed_dict:
                new_row[1] = passed_dict[obj_id]
                cursor_update.updateRow(new_row)
def main():
    update_points(build_dictionary())

main()


texts: Creating regular point layer within several buffers using QGIS?

texts:
I have several buffers distributed along a whole country (therefore, the distances between them is very high). These buffers are represented in the same layer/shp.

I want to create a regular point layer (spaced between 100m) only within each buffer, but so far the only procedure I know is to create a regular point layer in all surface, and then use "clip" tool to get only the points within each buffer.

However, it seems impossible (at least with my computer) to create a regular point layer for a whole country separated by 100m, and then clip it.

Is there any other way to obtain the regular point layer only within each buffer?


texts: How to add custom basemaps to BasemapGallery

texts:
I would like to add my own custom basemaps to the BasemapGallery. Is there any way I can create my own BasemapGallery with only my maps? I am using the ArcGIS API for javascript.

EDIT: for clarity

I have already created a custom basemap that the answer here describes. The issue is that I would like to create a Basemap Gallery but with my own basemap, not ArcGIS's. I hope this makes a little more sense now.


texts: I have  set a bookmark within QGIS and would like to set another bookmark just slightly east of the current one but with the same scale and north-south extent. If I use the pan tool the north-south extent changes because I can't hold my hand that still. I don't see scroll bars in the window and the arrows on the keyboad go too far to the east. Is there any other way I can shift the window slightly to the east without changing the latitude?


texts: Exact method used in simplify command within "Convert ArcGIS Raster to Polygons"

texts:
I am trying to understand the exact method used in the "simplify" command under the "Convert Raster to Polygons" in ArcGIS.

The help page describes qualitatively what the command does, but does not describe the method in detail. I have looked at other help pages on "simplify" polygons and come across a number of methods that can be used when simplifying polygon shapes, but it is unclear which of these methods are used in the exact Raster to Polygon command.

One that I found helpful was http://desktop.arcgis.com/en/arcmap/10.3/manage-data/editing-existing-features/about-simplifying-and-smoothing-a-feature.htm, which suggested to me that perhaps the "simplify" command in Raster to Polygon uses the "Generalize" that uses the published algorithm (Douglas and Peucker, 1973).

Is this correct?


texts: Google Earth Engine memory capacity exceeded

texts:
I am trying to export to a csv table with total precipitation per month-year per polygon (~5571 polygons) from 2001-2019 (228 time steps). I am using the CHIRPS dataset (daily data) so had to reduce the data to mean month-year.
I can download about two months of data but when I try to download 2+ years of data I get the message: Error: User memory limit exceeded.
I have tried every suggestion I have found searching online.
Is there a way to rewrite my code (linked here) to not get this error message?
// Set years and month
var startYear = 2001;
var endYear = 2019; //need to run in subsets?
var years = ee.List.sequence(startYear, endYear);
var months = ee.List.sequence(1,12);

// load the image collection
var Daily = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY")

//~5571 features (polygons)
var municipalities = ee.FeatureCollection("shapefiles")

//rescale features to help with memory error
var municipScaled = municipalities.map(function(feature) {
  return feature.simplify(100);
});

// make monthly summed mosaics
// loop over the years and months to get summed monthly images
var byMonth = ee.ImageCollection(ee.FeatureCollection(years.map(function(y){
  var yearCollection = Daily.filter(ee.Filter.calendarRange(y, y, 'year'));
  var byYear = ee.ImageCollection.fromImages(
    months.map(function(m) {
      var summedImage = yearCollection.filter(ee.Filter.calendarRange(m, m, 'month'))
                  .reduce(ee.Reducer.sum());
      var date = ee.Date.fromYMD(y, m, 1).format("MM_dd_YYYY");
      return summedImage.set('system:time_start', ee.Date.fromYMD(y, m, 1)).rename(ee.String("summed_precip")
                                       .cat(date));
      //.set('month', m).set('year', y); // eventually set year and month
  }));
  return byYear;
})).flatten());
print(byMonth)

// filter the empty one out
var outputMonthly = byMonth.filter(ee.Filter.listContains('system:band_names', 'constant').not())
                    .sort('system:time_start').toBands();
print(outputMonthly);

//test to make sure range makes sense
//Map.addLayer(outputMonthly.select("0_0_summed_precip01_01_2001"), {min:0, max:50});

//determine scale
var scale = Daily.first().projection().nominalScale();
print(scale)

//reduce to total precipitation per municipality
var muncip_monthly_precip = outputMonthly.reduceRegions(municipScaled, ee.Reducer.sum(), 5000, 'EPSG:4326')
                                .map(function(feature){
                      return(ee.Feature(feature).setGeometry(null)); // map over the feature collection and drop the geometry for memory saving
                    }).copyProperties(municipScaled, ee.List(["CD_MUN"]));

// save the table to google drive
Export.table.toDrive({
  collection: muncip_monthly_precip,
  description: "total_monthly_precip",
  folder: 'VL_GEE',
  fileFormat: 'CSV'})


texts:
There is plenty of remote sensing datasets out there, but I haven't been able to find one along with field data. Particularly I am looking for airborne hyperspectral images that have field data identifying plant species as well as their coordinates in the image.

Even if a dataset is from one organization and image from another, that works too.


texts: Postgresql Query Plugins Not Returning Results in QGIS

texts:
I'm using the Fishnet function as described here to create a grid.

If I use the following SQL to call the function to position a grid where I want it (I'm using NZGD2000 SRID = 2193) the query returns from postgresql as expected:

create sequence temp;

select
    nextval('temp')::int as id,
  st_setsrid(geom,2193)
from
    st_createfishnet(100,100,1000,1000,1645375,5349447);

If I try to create a layer in qgis using the select portion of the above SQL statement (i.e. not the create sequence command), using either 'Fast SQL Layer' or 'RT SQL Layer' plugins, the query returns zero results.

If I execute the following SQL statement in postgresql to create the temporary table and then add a PostGIS layer to qgis, the grid shows fine with all records returned as expected.

select
  nextval('temp')::int as id,
  st_setsrid(geom,2193) as geom
into
    temptable
from
    st_createfishnet(100,100,1000,1000,1645375,5349447);

I've used both of these plugins before for doing on-the-fly calculations in postgis and displaying in qgis - usefully eliminating the need for temporary tables. Is there a bug in them that is being triggered by this specific query? Or am I doing something wrong that I just can't see?

Using Quantum GIS 1.8.0 ; Postgresql 9.1 ; PostGIS 2.0

EDIT: clarified that I'm only running one statement in the plugin, without the create sequence command


texts: Reading text-anchor and text-angle from DXF using geopandas -> fiona -> gdal

texts:
I am working on a Pyhton desktop app that publishes CAD-files (dxf) to Postgresql/ GeoServer. It uses geopandas to import the dxf-file. The resulting geodataframe doesn't contain the anchorpoint and angle of the texts in the dxf-file, which I need to correctly display the texts in the GeoServer layers. I find very little information about this subject on the internet.
Geopandas uses Fiona to import the data and Fiona uses Gdal to import a dxf-file.
Does anybody know a way to read text-angles and text-anchorpoints from a dxf file using geopandas?
As an alternative route I am looking at using the python lib 'ezdxf'. I'm sure it's possible to solve the problem using this library, but if it is possible using geopandas the solution is likely to be much simpler.
Any tips?


texts: How does Google S2 / Hilbert Curve handle the outliers?

texts:
Regarding Google S2 library, it used for a quick indexing on locations.
E.g. near by locations would have similar indices. Thus it's much easier to make location queries.

But, as we map a 2-D point to 1-D point, it's bound to lose some precision.
E.g. the following case, points are close but indices are far away.

Is there a pattern for the outliers and how does google s2 cope with it?


texts: Weight from field option not showing all fields from attribute table in heatmap algorithm in QGIS 3.18.1

texts:
I'm trying to create a heatmap on a shapefile of numerous points based on a single field in their attribute table.
The header of the attribute table is this:

... and I want to make a heatmap weighted on the "frp" field from the attribute table. However, in the Heatmap (Kernel Density Estimation) settings, I cannot see the "frp" field - only the "longitude" and "latitude" values appear:

What could be wrong here?


texts: QGIS plugin - NameError: name 'iface' is not defined while importing class from other file

texts:
Here is my class

class aClassName():
    def __init__(self):
        #do something...

        iface.setActiveLayer(QgsProject.instance().mapLayersByName('someName')[0])

        #do something
    ...

The class calls a number of methods from iface such as setActiveLayer, .mapCanvas().zoomToSelected() etc in different functions. It runs properly when it is called in the same file.

But when I tried to use the class by importing it into another file:

from theClassFile import aClassName

The error message showed up and said :

  NameError: name 'iface' is not defined

even if adding from qgis.utils import iface to theClassFile with reference from here.

I just began to use PyQGIS last month.

Here is my QGIS version:


texts:
I´m trying to build a polygon from a list of coordinates. I have the coordinates in different ways as:

pointsfloat
[(583899.8, 8177297.19), (583900.23, 8177262.52), (583896.49, 8177262.47), (583895.96, 8177297.29), (583899.8, 8177297.19)]

or

coordenadas
[583899.8, 8177297.19, 583900.23, 8177262.52, 583896.49, 8177262.47, 583895.96, 8177297.29, 583899.8, 8177297.19]

When I try to convert it to a QgsPoint I always lose the decimals digits, it's converting the numbers to integer automatically.

How can I generate the polygon feature from those vertices? I'm using python 2.7 and qgis 2.14.


texts: Does merging segments in Google Map Maker ever work?

texts:
When editing line segments in Google Map Maker, the number of crossings one needs to touch for a change is pretty relevant. From my experience, every move of a crossing needs to be reviewed, so the number of crossings needed to be moved is the number of round trips one needs for a change.

So, one way to improve this would be to remove crossings which anyway only connect two segments first. However, I've never been able to do this. I always get an error message saying that the segments to be merged have different attributes. I tried to find the differences, but without success.

But let's start with a general question: Did anyone ever successfully remove a crossing in Google Map Maker?

This is the full error message (in German):

  Segmente mit ungleichen Attributen können nicht zusammengeführt und
  die entsprechenden Kreuzungen können nicht gelöscht werden. Die
  Segmente, die an dieser Kreuzung enden, weisen Attribute auf, die
  nicht zueinander passen. Bitte korrigieren Sie dies und versuchen Sie
  dann erneut, die Kreuzung zu löschen.

and it translates to:

  Segments with unequal attributes can not be merged and the
  corresponding intersections can not be deleted. The segments ending at
  this intersection have attributes that do not match. Please correct
  this and then try to clear the intersection again.


texts: Mapshaper CLI: Trying to filter on iso_a2

texts:
I'm new to Mapshaper and I tried to find a good tutorial, but found none that answered my questions.

I'm using mapshaper on the commandline
using Natural Earth Data
want to use filter to select only the a specific group of countries

There is a page with an example of selecting only a few counties from a state, this I adapted to:
mapshaper ne_10m_admin_0_countries_lakes.shp -filter '"NL,BE,CH".indexOf(iso_a2) > -1' -o format=svg out.svg
I get:

[filter] ReferenceError in expression [iso_a2==FR]: FR is not defined

Running the following, simplified command:
mapshaper ne_10m_admin_0_countries_lakes.shp -filter 'iso_a2=="FR"' -o format=svg out.svg
I still get:

[filter] ReferenceError in expression [iso_a2==FR]: FR is not defined

Anyone out there with the tip to get me started in the Mapshaper game?


texts: How to query elevation data from "Digitales Geländemodell (DGM) NRW", maybe using WCS?

texts:
I use street network data originating from OpenStreetMap, stored in PostGIS, and want to add corresponding height/elevation data from "Digitales Geländemodell (DGM)", available for North Rhine-Westphalia/Germany from "Geobasis NRW":
https://www.bezreg-koeln.nrw.de/brk_internet/geobasis/hoehenmodelle/gelaendemodell/index.html
This page lists several links to get data. As far as I understand only the downloadable xyz files and the first WCS link ("wcs_nw_dgm") contain the actual height data in m/cm.
"xyz" are huge text files, 124 MB for each square of 2 x 2 km, in total about 9000 files for North Rhine-Westphalia. I downloaded some of them and wrote a Python script to query elevation of highway nodes in them. This is working, but I would prefer to use less of download and disk capacity.
I found I can bind the WCS service in QGIS and can query single point elevations with QGIS's Info tool, so I hope this could be used to selectively query a few thousands of network nodes? But it seems WCS itself is not usable for point queries, but returns data for rectangles as tif/png/jpg image.


texts:
I have two shape files and want to join these two attribute table with common field. But arcpy throws the following error.

arcpy.management.AddJoin(output_feature, "AccountPoi", lstfc1[x1-1], "AccountID", "KEEP_COMMON")

Traceback (most recent call last):
  File "C:/Users/dxc163530/PycharmProjects/start_charm.py", line 48, in <module>
    arcpy.management.AddJoin(output_feature, "AccountPoi", lstfc1[x1-1], "AccountID", "KEEP_COMMON")
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\management.py", line 5096, in AddJoin
    raise e
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\management.py", line 5093, in AddJoin
    retval = convertArcObjectToPythonObject(gp.AddJoin_management(*gp_fixargs((in_layer_or_view, in_field, join_table, join_field, join_type), True)))
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 506, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
The value cannot be a feature class
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer.
WARNING 000970: The join field AccountPoi in the join table AccountPoints is not indexed. To improve performance, we recommend that an index be created for the join field in the join table.
Failed to execute (AddJoin).

Process finished with exit code 1


texts: keeping the coordinate system of raster files in the resulting raster file after operation with numpy

texts: I'm trying to read as array two raster files using ReadAsArray in GDAL/python, perform some map algebra operation on it using numpy and then write the resulting raster file in GTiff format using WriteArray in GDAL/Python. The spatial reference system of the input raster files are the same. I want to know the resulting GTiff file will have the same coordinate system or not?


texts:
I have a bunch of thousands of points.
I'd like to group them, in a user-defined number of clusters for example or based on some kind of "best number of clusters", based for example on their relatives distances.
Which kind of algorithm is available on postGIS and/or shapely (preferred) to do that?
I can simply make some buffers around them and a union but I wonder if there is some other advanced technique that are computationally efficient on large dataset.

Edit:

I've found the DBSCAN algorithm [1] interesting and it seems rather simple to insert in an existing code. There is also the OPTICS algorithm [2] but it seems harder to find a solution.

[1]https://en.wikipedia.org/wiki/DBSCAN
[2]https://en.wikipedia.org/wiki/OPTICS_algorithm

Anyway, I'm stuck with DBSCAN now:
how to store the cluster label in which a point lays within a new column "cluster" in a postgreSQL table? Actually, the points ID are not stored along their coordinates within the cluster...
Some checking operations on coordinates in order to retrieve points IDs in original array would be a bit overrated, in my opinion.


texts: PTV_Mercator correctly recalculated to WGS84 EPGS:4326 but not correct on map

texts:
I have recalculated PTV_Mercator (PTV Spherical Mercator EPSG:505456) coordinates to WGS84 EPGS:4326 coordinates in Excel. And when checking them on Google Maps they show up correct.

But after loading the CSV and trying to combine them with a map (tried several like OSM, QGIS, etc.), I can in no way get them to "overlap". I've read (and tried) the QGIS "Projection"-topic in the manual several times, but I can not get it to work correct.

EDIT:

As requested, I have added a picture of the problem:

The recalculation was done with the use of this code-snippet:


texts: Removing space between multipart polygon in QGIS 3?

texts:
I am looking for an efficient method to edit and remove the space between multipart polygons in QGIS 3. I have an example below which should hopefully be clear:

Multipart polygon

Final result

I know I can trace around the perimeter to create a new feature but this isn't always practical. Ideally I would like to edit it in place without the need of creating a new layer.


texts:
I'm trying to create a bounding box from a simple Lat/Long, given an increase size in metres.

For example, create a bounding box from X/Y which goes outwards 5000m at the NE / SW corner. (yes, this isn't a circle so the radius isn't even, etc).

Result:

NorthEast corner Lat/Long
SouthWest corner Lat/Long.

My guess is I need to convert meters to ..something? are they Radians?

I'm used to just using SQL Server geospatial to create new polygons from a point. e.g. @somePoint.STBuffer(5000) creates a new circle with a radius of 5000 meters, from the centre point (given the point is WSG84).

How can I do this?

EDIT 1:

As requested re: GIS tools - this is tools agnostic. As in, i'm just after some common math formula to calculate the new co-ordinates. I'm not using any GIS program, but doing some coding.

EDIT 2:
- This will be used to search for all points in a poly, on Earth. I just need to create a rough/simple square from a Lat/Long (GPS coordinate).
- Then finally, I will either do this in a spatial query in a DB (for example, RavenDb, SqlServer, PostGIS, etc) or more likely, in .NET C#. So I was hoping to see if there was a common algorithm  before I checked if this has been coded up in some coding language.


texts: Optional Input and Workspace in ModelBuilder?

texts:
1 - I am trying to create a model to buffer rivers, lake and spring but sometimes I don't have some of my input (river, lake or spring). Is it possible to create a model with optional input ? How to do that ? (I don't have experience with Python).

Another question is:
2 - We work whith 5 computers and I want to share the model with others (We don't have server, every user has its own user name) How should I set the workspace and how to do that ?


texts:
I have a couple of shapefiles, the first one is a set of points, while the second one is a set of lines. I need to regularize the shapefiles in order to be able to generate a graph (adjacency list / paths) out of them later.

The set of lines has many dangles, with both undershoot and overshoot problems. I want ArcGIS to mark the suspect points and possibly offer me a correction to apply, on an individual basis.

I've seen the tutorial "Using geodatabase topology to fix line errors", which seems to go where I want.

However, in order to create a "topology" with a "Must Not Have Dangles" rule, I need a "geodatabase".
I'm very new to GIS systems, and I'm trying to figure out a way to convert my shapefile to a geodatabase. I did the following

open ArcCatalog
click Connect Folder to connect the folder with my shapefiles
right click on connected folder, New, File Geodatabase
right click on created geodatabase, New, Feature Dataset
right click on created feature dataset, New, Feature Class
added the shapefile with lines (I used the Import button while creating the new Feature Class - see first screenshot below).
right click on the feature dataset (not the feature class), New, Topology
follow the tutorial

However, after adding the "Must Not Have Dangles" rule, and running Validate, I am told there are no errors. I can see there are errors, there are obvious dangling lines, so why aren't they detected?

When I'm finished creating the topology, I get this warning:

This is the final result:


texts: Extract by mask: area difference in output result?

texts:
When we do extract by mask in ArcMap results always shows some area difference in raster outputs.

We are using clip layer as vector and layer to be clipped as raster DEM.

Our area of interest is irregular vector boundary but when we extract those areas from raster DEM, the output(extracted DEM) area will not matching with actual vector area we used.

How can we solve this problem?


texts:
I'm a novice with Carto and GIS generally.  I have a data table with a few thousand lines of data.  Each row of data has a country and a numerical value.

I'm trying to total the value by country and show a single point per country.

I've tried reading documentation, but am overwhelmed by reading up on PostgreSQL - what I am attempting to do is so simple, yet I don't find a way of making it work.

Should I change my data structure - to only have one line of data per country?  I would be surprised is that was necessary!

As you see in the map - the Netherlands has 3 or 4 overlapping circles - there should only be a single circle.

Edit: additional picture.


texts:
I am trying to find a tool in QGIS that does the same function as Region group in ArcGIS.

What I am trying to do is to identify uniquely each patch of a raster that has a certain value (e.g "Value" =1) using neighbouring cells.

From the region group docs:

  For each cell in the output, the identity of the connected region to which that cell belongs is recorded. A unique number is assigned to each region.


texts: Obtaining time series of water extents per plot using landsat NDWI data available on GEE

texts:
I wanted to obtain a monthly time series estimate of the water surface for each of 67 survey plots (500m circles). However, I could not figure out how to correctly code this in the google earth engine. My code below does not seem to count the number of water pixels per plot.
Here is my code link: https://code.earthengine.google.com/b62cd66b38e974cfbb402100dc10e3bd
// // Load from Fusion table
var surveyPlots = ee.FeatureCollection(table, 'geometry');

// Create water mask
var water = function(image) {
  return image.updateMask(image.gt(0));
};

// Load a dataset.
var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_8DAY_NDWI')
  .filterBounds(surveyPlots)
  .filterDate('2014-01-01','2020-08-01')
  .map(water)
  .select('NDWI');

var bandName = ee.Image(dataset.first()).bandNames().get(0);

var startDate = ee.Date('2014-04-01'); // set analysis start time
var endDate = ee.Date('2014-10-01'); // set analysis end time

// calculate the number of months to process
var nMonths = ee.Number(endDate.difference(startDate,'month')).round();

// get a list of time strings to pass into a dictionary later on
var monList = ee.List(ee.List.sequence(0,nMonths).map(function (n){
  return startDate.advance(n,'month').format('YYYMMdd');
}))
print(monList)

var result = surveyPlots.map(function(feature){
  // map over each month
  var timeSeries = ee.List.sequence(0,nMonths).map(function (n){
    // calculate the offset from startDate
    var ini = startDate.advance(n,'month');
    // advance just one month
    var end = ini.advance(1,'month');
    // filter and reduce
    var data = dataset.filterDate(ini,end).median().reduceRegion({
      reducer: ee.Reducer.count(),
      geometry: feature.geometry(),
      scale: 30
    });
    // get the value and check that it has data
    var val = ee.Number(data.get(bandName));
    val = ee.Number(ee.Algorithms.If(val,val,-999));
    // return value with a time key
    return val;
  });
// create new dictionary with date strings and values
  var timeDict = ee.Dictionary.fromLists(monList,ee.List(timeSeries));
  // return feature with a timeseries property and results
  return feature.set(timeDict);
});

// print to see if it is doing what we expect...
print(result);

// drop .geo column
var polyOut = result.select(['.*'],null,false);

// Export the data to a table for further analysis
Export.table.toDrive({
  collection:polyOut,
  description:"time_series_NDWI",
  fileFormat:"CSV",
});


texts: Geotagged photographs hosted on website and shared via KML

texts: I have a database of about 500 geotagged photographs that I would like to host on my personal website, and then share as a KML to other users. They would get a relatively small KML file size (~1-2 MB) with placemarkers of each photo and thumbnails. When you click the thumbnail the path would direct it to the webpage to pull up the full photo in a browser. I have access to Google Earth, Arcmap 10.2, Matlab, etc. Seems Python would be the ideal tool for this but I don't have much experience there. Can anyone point me to the general steps I would need to make this happen?


texts: Point (of a linestring) within Polygon using ogr and Python

texts:
I'm currently working on a project in which I need to build a topological network out of the geometry features I find in shapefiles.
So far using Ben Reilly's open source project I've managed to transform linestrings into networkx edges, as well as detecting close features (other linestrings say) and adding them to the closest point so that I can run shortest path algorithms.

But that's fine for one shapefile. However, I now need to connect features from different shapefiles into a big networkx graph. So for instance, if a point is within a polygon I would connect it (by connect it I mean add a networkx edge - add_edge(g.GetPoint(1), g.GetPoint(2)) with the point in the next shapefile that also is within a polygon that shares a similar attribute (say, ID).
Note that the polygons in the different shps only share the same IDs and not the coordinates. The points that fall within the polygons also do not share the same coordinates.

My solution to this problem was to identify the point that resides in a polygon, store it, find the point in the next shapefile that resides in the polygon with the same id and then add the networkx edge between them.

How to find if a point resides within a polygon? Well, there's a well known algorithm: RayCasting algorith that does that.
This is where I actually got stuck though, because in order to implement the algorithm I need the coordinates of the polygon, and don't know how to access them right now, even after skimming through a documentation of the OGR's Geometry.
So, the question that I'm asking is how to access the polygon points, or coordinates OR is there an easier way of detecting whether a point falls within a polygon?
Using python with osgeo.ogr library I coded the following:

 if g.GetGeometryType() == 3: #polygon
                c = g.GetDimension()
                x = g.GetPointCount()
                y = g.GetY()
                z = g.GetZ()

see the image for a better understanding of my issue.

[EDIT]
So far I've tried storing all the polygon objects in a list with which I would then compare the linestring first and last point. But Paolo's example is related to using the Point Object reference and the Polygon Object reference, which would not work with a line object reference since not the entire line is within the polygon but rather the first or the last point of its linestring.

[EDIT3]
Creating a new Geometry point object from the coordinates of the first and last point of the linestring and then using that one to compare against the polygon geometry objects saved in a list seems to work just fine:

for findex in xrange(lyr.GetFeatureCount()):
    f = lyr.GetFeature(findex)
    flddata = getfieldinfo(lyr,f,fields)
    g = f.geometry()
    if g.GetGeometryType() == 2:
        for j in xrange(g.GetPointCount()):
            if j == 0 or j == g.GetPointCount():
                point = ogr.Geometry(ogr.wkbPoint)
                point.AddPoint(g.Getx(j),g.GetY(j))
                if point.Within(Network.polygons[x][0].GetGeometryRef()):
    print g.GetPoint(j)

Thanks Paolo and Chris for the hints.


texts: Check if point is within radius of another point with Shapely

texts:
I've successfully used Shapely with python to check if a point is within one of a set of polygons in a .shp file but now would like to check if a point is within a circle of radius x miles from another point - how can I achieve this with Shapely?

I'd imagine I have to create a circle with radius x miles from my reference point and then check if each of my other points intersects with that new circle - I'm ok with checking the intersection but how do I create a circle with radius x miles given my coordinates are in lat/long?

I read a question that suggests I need to change projections but not sure what this means or how to apply it to my problem: creating circle with radius in metres

My example is in Africa btw, not the US (the above question mentions the US local state plane as a 'reference system' but not sure what that means)...


texts: Openlayers 3 equivalent of Openlayers 2 OpenLayers.Protocol.HTTP.options.params?

texts:
Is there a way to send data to the server with POST request instead of GET ?

I used to use Protocol.HTTP.options.params for this purpose, but I can't find any equivalent in Openlayers 3 OpenLayers.source.Vector.


texts:
Is it possible to get mapbasic to use the current directory for processes?

So rather than the following where WFolder is manually set as the working folder can I get it to just work in the folder where the mbx is?

Include "MapBasic.def"
Include "Menu.def"

Declare Sub Main()
Dim filename, printtext,WFolder as string
Dim rcount,i as Integer

Sub Main()
Close All
'OnError GoTo CaughtEx '// Turn Error handling on
WFolder="D:\Projects\AsCon_data\Output\OceanBreeze_4B_4C\"
Register Table "D:\Projects\AsCon_data\Output\OceanBreeze_4B_4C\List.txt"  TYPE ASCII Delimiter 9 Titles Charset "WindowsLatin1" Into WFolder+"~MAP0001.TAB"
Open Table WFolder+"~MAP0001.TAB" Hide
Commit Table ~MAP0001 As WFolder+"List.TAB"
Close Table ~MAP0001
Open Table WFolder+"List.TAB"  as List Interactive
Browse * From List
rcount=TableInfo(List,TAB_INFO_NROWS)

i=1
While i<=rcount
    Select * from List where List.RowID=i into RSelection
    Browse * from RSelection
    print RSelection.col1
    Open Table RSelection.col1 as temptable

' Do some action(s)

    commit table temptable
    i=i+1
    close table temptable

Wend
CaughtEx:  '// error handling routine
    '// Log error here...
    print Chr$(13)+Chr$(13)+"///////MAIN///////"+Chr$(13)+"FAILED due to:"+Chr$(13)+Error$()+Chr$(13)+"\\\\\\\\\\\\\\\"+Chr$(13)
    Resume Next  '// resume code at line after the one which threw the error
End Sub


texts: How to calculate overlap area of two polygons?

texts:
I am a beginner in the world of GIS.

I have one table polygons with each polygon having north-east and south-west coordinates. And for every polygon, I need to find area of overlap with other polygons (e.g. 70% area of polygon A overlaps with polygon B).

I know the ST_Intersects function, but it returns boolean. What other function that I can use with putting buffer as 70%?

If anybody could get me started, that would be great!


texts: WFS Layer not displaying in ArcGIS and QGIS, but appears in the browser

texts:
I am using django-gisserver to open GIS models to WFS layers. I have configured the view and URLs in my Django project without any issues, and when I check the WFS URL running our Django project locally, the WFS and the associated features work, and the geometry fields are populated. However, when I then go to QGIS and ArcGIS to test the functionality, the layers import but do not display. The attribute table is populated though.
Has anyone else experienced this setting up WFS layers?
edit: here is the URL I am using: http://127.0.0.1:8000/wfs/unitgeos/. I am working in a shared repository and have code review to pass before I can share a public-facing URL.


texts: OpenStreetMap how to identify major ports

texts:
I'm quite new to OSM data.  Can anyone help me with the best tags to use to identify major shipping ports?  I'd like to be able to identify cities that have a significant amount of commercial shipping activity.  I haven't found this as specific tags, so maybe I need to use something else as a proxy for shipping activity.  I don't know - something like  cities that are harbours and have customs and a lot of docks and cranes, etc etc.

I'd love something a little simpler, if you're able to help.


texts: return when asynchronous task has completed

texts:
I'd need a function to wait until the asynchronous task has completed and then return.

In my example, I'm getting undefined returned from the "execute" method:

var result = myClass.execute(obj);  //undefined

How can I get the actual result returned?

return declare(null, {
        constructor: function (obj) {
            this.mapServiceUrl = obj.mapServiceUrl;
            this.geometry_json = obj.geometry_json || undefined;
            this.where = obj.where || "1=1";
        },

        execute: function () {
            var queryTask = new QueryTask(this.mapServiceUrl),
                query = new Query()
            ;

            query.returnGeometry = true;
            query.outFields = ["RDUWI"];

            if (this.geometry_json) {
                query.spatialRelationship = Query.SPATIAL_REL_INTERSECTS;
                query.geometry = Polygon(this.geometry_json);
            };

            query.where = this.where;
            var deferred = queryTask.execute(query);
            deferred.then(function (featureSet) {
                return featureSet;
            });
        }
    });


texts:
I started recently working in a project where I need to start using EPSG codes to identify coordinate systems and their elements.

I saw that a single EPSG code can identify multiple objects. I am interpreting that each object must be from a different type (i.e., no two Datums can have the same EPSG code, but a Datum and a CRS can).

Am I correct with that assumption?


texts: Input of in-memory geometry object to Dissolve tool in ArcPy script throws Error 000732?

texts:
In an arcpy script, I run a Union on a list of geometry objects.  The output is set to an empty geometry object.  The Union runs fine.  What I am trying to do is then feed that output into the Dissolve tool.  No matter how I format the input, I get the following error:

Parameters are not valid.
ERROR 000732: Input Features: Dataset in_memory\f8537CACB_AB51_4BDC_B8E7_0895F11E6AF0 does not exist or is not supported

ESRI help files indicate that tools set to output to an empty geometry object return a list of geometry objects.  The output from the Union tool gives <Geometry object at 0x26d25910[0x27397068]> when queried with repr().  I have tried putting the Union result directly into the Dissolve tool and likewise the 0 index item from the result, seeing as it supposed to be a list.  Both result in an error, the latter saying that a geometry object is not indexable.  Here is a code snippet:

if len(lstHSGeoms) > 1:
            #Union geometries into one geometry
            lstHSUnionGeom = arcpy.Geometry()
            arcpy.Union_analysis(lstHSGeoms, lstHSUnionGeom, "ONLY_FID")
            arcpy.AddMessage('          lstHSUnionGeom:  %s' % repr(lstHSUnionGeom))
            lstHSOneGeom = arcpy.Geometry()
            arcpy.Dissolve_management(lstHSUnionGeom,lstHSOneGeom, None, None, "SINGLE_PART")
            hsOneGeom = lstHSOneGeom[0]

The list lstHSGeoms is a Python list of polygon geometries, which have checked out as valid.  The Union runs, but the Dissovle doesn't.  I'm stymied.  BTW, I tried using an empty string for the fourth parameter of Dissolve (seeing that it is supposed to be a string type), but it didn't make a difference.


texts: How can I get the default zoom level to work in QGIS for geography types?

texts:
Create a table like this, and then add it as a layer to PostGIS. It has one simple geography, a line.

CREATE TABLE foo AS
SELECT ST_GeogFromText(
  'LINESTRING(-122.27172 37.80345,-122.27172 37.81549)'
) AS highway;

You can even create an index on it if you think it matters,

CREATE INDEX hwidx ON foo USING gist ( highway );

You can see the line in WKT in the table,

                      st_astext
-----------------------------------------------------
 LINESTRING(-122.27172 37.80345,-122.27172 37.81549)

Here is the table definition,

       Table "public.foo"
 Column  |   Type    | Modifiers
---------+-----------+-----------
 highway | geography |

Yet, when I Zoom to Layer in QGIS, no joy! It just shows a blank screen. I am using QGIS 2.18.2.


texts: Trying to extract city data of city lights location?

texts:
Here is a website:
http://geodesy.net/webmaps/milpitas/gedit2/gedit.htm?map=base

I believe they are using Google Map API and Fusion Table(deprecated 2016).

I just wanted to extract the data to a text format.


texts:
I just had a really hard time getting apt to add the repository key to a fresh install of Debian 10.

After much trial and error, I discovered that the following steps work:

wget https://qgis.org/downloads/qgis-2017.gpg.key
sudo apt-key add qgis-2017.gpg.key

When I tried the commands listed on the install page, I got this error:

gpg: key 51F523511C7028C3: 1 signature not checked due to a missing key

Also, note that I only had success with the 2017 key; not the 2019 one.


texts: Using feature service with JavaScript API and ArcGIS for Server on localhost?

texts:
I installed ArcGIS 10.2 for Server on my localhost and created and published a feature service on it. Now I want to use the feature service I created on a simple html page using Javascript. How can I do that?

I read the resources on arcgis.com but it just described how to do it with the maps on arcgis.com!

I don't need a base map from arcgis.com. actually I want to use the web page offline(just on localhost).


texts: Join (map-matching) observations to roads

texts:
I know that the subject is wide, but I could not find answers to my problem. I have to shapefiles: Roads(c.a. 580 rows and divided into three categories of roads) and observations (detected speed limit signs - c.a 11000 rows). I have to join (map-matching) observations to roads.
I've already tried Spatial Join (match option - closest) in ArcGIS Pro and NNJoin plugin in QGIS, and both gave me some results, but I am still not happy.
The first issue is that the layers have only a geographic coordinate system, and I wonder if they should also have a Projected coordinate system to perform this join?
The second issue is the fact that when using the closest approach some of the points can be assigned to the wrong segment of the road which can give me some errors in the output. I cannot use any other match option because the two layers do not intersect with each other, whatsoever.
How could I approach this issue differently? I've already checked this post:  Map-matching. Join points to lines, but this not fully answer my questions.