texts: Converting xy coordinates in ITRF2000 to long lat WGS84

texts: How can I convert the x y coordinates which are in ITRF2000 coordinate system to long lat WGS84?


texts: Automating entire map production in ArcMap using ArcPy/ArcObjects?

texts:
I am trying to gather some information on how to create maps based on an existing base map I have. The basemap has all the data loaded. So the user will enter data into a form and the  would do a save as, query data in the map, create the map and export to PDF.

I know python is for automation of workflows, but what I need is a way to have my users enter information into a form and have the map created automatically based on a save as of an existing map, some inputs, queries and data in my basemap. I see lots of examples for automating workflows, but I am looking to deploy a solution that non-GIS users could create maps quickly based on some information in arcmap.

Would I do this with ArcObjects or can it in fact be done with Python? Has anyone tried automating the entire map process?

I just need to know which method I would use to create this project? And if I can complete with Python and/or ArcObjects?


texts: Deciphering coord. ref. in R – very basic question

texts:
I'm trying to make sense of the coord. ref. parameters for rasters in R. For example, for a given raster, crs() gives the following:

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

The datum and ellps parameters make sense, but the proj argument has me confused. Does this mean that the raster is unprojected as opposed to something like +proj=utm +zone=36 which would refer to projected data? Also, what is the meaning of +no_defs? Unfortunately, I have been unable to find documentation that explains the string output by crs().


texts: How to display the Radial type of pattern associated with electric consumers

texts:
I have a point dataset of electrical poles (feature data set under FGDB) and also i have the consumers data in non-spatial format (excel sheet), the excel sheet having a pole ID field which is the common  field with point dataset as well. I have joined the data and generated an output file, which having all the data in it.
example - Pole ID -ABCD/1 showing 5 Consumers joined with it.
Here i want to display the radial type of pattern of consumers connected with the pole only when i click that particular pole (Point). I couldn't find any such solution for in ArcGIS.
Below attached image will illustrate how i want to show the data.
Example data attached...


texts: ArcGIS API for JavaScript - check if a polygon exists in a polygon feature layer

texts:
I have a polygon feature layer and a polygon geometry. How do I find out if the feature layer contains the polygon, as in exact shape of the polygon?

I did see the geometry engine, but that doesn't contain a method for knowing if exact shape exists at same place. Both datasets have same spatial reference.


texts:
I am designing a map of very densely populated area in Tilemill. I am using OSM building data. Due to having a lot of buildings it is very difficult to differentiate one building to another.

Is there any way to reduce all building polygon sizes to 95% or 98% so that buildings can easily be differentiate each other?


texts:
I have this function works totally fine. But i want to be able to change the font size

 def label(wah, field):
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.name == wah:
            if lyr.supports("LABELCLASSES"):
                for lblClass in lyr.labelClasses:
                    lblClass.expression = "["+field+"]"
                    lyr.showLabels = True

I have tried something like this but i cannot quite seem to get it to work. The label size is not changing from the code below

 lblClass.expression = "<FNT size = '8'>" + "[" + {} + "]" + "</FNT>".format(field)

or

lblClass.expression = '"{}" + [{}] +  "{}"'.format("<FNT size = '8'>",field,"</FNT>")

Just for reference when i am calling the function it looks like this:

label("CountyClip", "County")

Can someone help me fix my code?


texts:
I am trying to understand how to create a database view that shows a 150 foot buffer around each point in a table.  The table of points changes pretty frequently and I would like to use a view instead of just buffering the points and having to update the polygon layer all the time.

I am using ArcGIS 10.2.1 and my SDE database is running on SQL Server 2012.

I have some experience building view layers in a PostgreSQL environment but that was a while ago and I'm struggling to get my head around doing things in this new space.

I started out reading the documentation for the 'Create Database View' tool but am having trouble figuring out how to put together the SQL statement.  Where do I put the .ST_Buffer() part of the statement?  I'm confused because .ST_Buffer() only accepts the distance as an argument, and it has a dot prefix, which makes me think that I need to use it as a function. I can't seem to find any references to doing this anywhere on the web.  Is it even possible to do this in the software environment I'm using?

This is the best stab I have taken:

SELECT OBJECTID, LocationID, shape.st_buffer(47.5) FROM gisdb.DCGIS.TransformerBank

When I run that the tool errors out and reports:

ERROR 999999: Error executing function.
[[Microsoft][SQL Server Native Client 11.0][SQL Server]Cannot find either column "shape" or the user-defined function or aggregate "shape.st_buffer", or the name is ambiguous.]
DBMS table not found
Failed to execute (CreateDatabaseView).


texts: How to split a line into a set of equidistant points

texts:
Used to rely on XTools for this, but in my current environment, I don't have access to it.

I'm trying to do this in ArcMap10.

The goal is to break a line representation of a stream into a set of equally-spaced points in order to eventually determine z values at those points so that I can determine gradient.

My intended workflow is as such.

split the stream into points
intersect the points with a DEM using GME
By determining the distance between the points, I should be able to determine that the gradient is.
Seems like it would make sense to integrate the point data back into the original streamlines. Unsure about this step.

I appreciate any critiques of my method, but my priority at this point is converting those streams to points.

Thank you!


texts: OpenLayers - vector layer not redrawing correctly when map zoom

texts:
I am loading an external GeoJSON file into a vector layer as described in the doc.

I load it into the map normally:

geojson_layer = new OpenLayers.Layer.Vector("GeoJSON", {
        strategies: [new OpenLayers.Strategy.Fixed()],
        protocol: new OpenLayers.Protocol.HTTP({
            url: "json/CS_Fr_St_General.json",
            format: new OpenLayers.Format.GeoJSON()
        })
    });

olmap.addLayer(geojson_layer);

When I zoom in or out, the old drawings of the geojson are not removed and the new one is added over it.

I suppose there is a setting in olmap that prevent the geojson to be refreshed/redrawn when zooming but I dont know how to fix it.

You can inspect the olmap and the geojson_layer in the debug tool from this address

https://carto48dev.mels.gouv.qc.ca/Infogeo/test/window.html

Use console.log(olmap) or console.log(geojson_layer)


texts: Selecting Matching Features in shapefile

texts:
I have a shapefile that draws up all 200 districts in a country. I want to select and colour 100 of these that I have on my list.

What would be the simplest way of going about this rather than me going through and selecting each one individually?


texts: Is it possible to inherit a session in ArcGIS 10.0 through the Python window?

texts:
So, I'm still struggling with an ESRI 10.0 installation with an Oracle ArcSDE database. I'm nearing the end ....

In trying to programmatically edit the versioned data. It needs to be done in an editing session. The DA module was not introduced to 10.1 so is unavailable. We thought maybe we could do it through the Python window in the Desktop environment. I have been trying but I don't believe there is a way to tell the Python script, use the current settings of the Desktop. I've tried some different combinations to get the ListFeatureClasses() but it is always empty.

I don't think I am going to be able to create / update / delete objects in this configuration. Or am I missing something?


texts:
In QGIS 2, when you click on a group which contains subgroups, all subgroups are checked/unchecked at the same time.
In QGIS 3, when you click on the same group, only that group is checked. Its subgroups remain unchanged.
I know there's an option in the submenu when you right-click the group, you can check/uncheck all of its children:

But is there a way this happens when you just click on the group instead of going through the menu? Is there a UI setting somewhere or some basic python code to achieve the same behaviour as in QGIS 2?


texts: Advantages of creating Python Toolbox over Custom Toolbox?

texts:
I enjoy making different geoprocessing tools to practice Python and to become more familiar with the arcpy modules. I've just been using a Custom Toolbox and then importing the Python script and setting the parameters, but I'm wondering if I should be using a Python Toolbox instead.

Are there any distinct advantages of one Toolbox over the other?


texts: Getting geometry properties of an ee.Image?

texts:
I try to extract the centroid of an ee.Image, as follows:

ee.Initialize()
image = ee.Image('LANDSAT/LC08/C01/T1/LC08_044034_20140318')\
          .geometry()\
          .centroid()\
          .getInfo()
# {'coordinates': [-122.14451407746256, 37.47160285059076], 'type': 'Point'} #Ok!

However, when I apply the same idea in the next example, I get:

ee.Initialize()
collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')\
               .filter(ee.Filter.eq('WRS_PATH', 44))\
               .filter(ee.Filter.eq('WRS_ROW', 34))\
               .filterDate('2014-01-01', '2015-01-01')\
               .sort('CLOUD_COVER')
image = collection.median()
image.geometry().getInfo()
#{'coordinates': [[[-180.0, -90.0], [180.0, -90.0], [180.0, 90.0], [-180.0, 90.0], [-180.0, -90.0]]], 'type': 'Polygon'} #Bad

I was expecting to get the bounding for the WRS_PATH and WRS_PATH specified, not the world bounding. Does someone have an idea to achieve this?


texts: In qgis cloud, how can I set what field is shown for attribute pop up with multiple layers?

texts:
I have a web map in QGIS cloud with several layers that intersect.  When I click on an area, I get the gray attribute pop up window showing me the layers that have features overlapping there, like this:

However, I want to be able to set the field the pop-up window shows. In the example here it showed "428" for Subdivisions which is the qc_id, but I would rather it showed "subdivisions". Under Parcels, it shows owner names whereas I would rather it show me Folio numbers.

Is there any way to change those fields to show different initial values?


texts: Drawing boundary of set of points using QGIS?

texts:
I have a set of points which I would like to change into a polygon.

I tried concave hull but do not get the intended polygon as you can see below,

And when I use Denaulay Triangulation I need to delete the outer triangles.

Is there an efficient way of doing the same?

I also tried Convex hull.


texts: Find minimum downstream flow distance from each polygon to any polygon in other feature class

texts:
The Problem:

I have two polygons feature classes, let's call them "SourceAreas" and "DangerAreas".  For each SourceArea, I need to know if water could flow to any DangerArea and if so, what is the minimum distance that the water would flow before it reaches any DangerArea.

Using ArcGIS Desktop 10.5.1.

What I Have Tried:

I have used a DEM to create a flow direction raster and then a flow accumulation raster with weights based on my SourceAreas (all other cells with a zero weight), and then a "Streams" feature class using the 'Stream to Feature' tool (Spatial Analyst -> Hydrology).

So now I can clearly see all possible water flows from all of my SourceAreas on the map, and visually I can see if they intersect any of my DangerAreas downstream.

However, there are a LOT of these and I need to automate the process, and to measure flow distances.  The attached screenshot illustrates a small sample of my overall area of interest.  The red areas are the SourceAreas, the yellow is the DangerAreas and the blue lines are the streams.

I have tried creating a geometric network from the streams feature class and was going to try to do a trace downstream from all the points where the streams intersect my SourceAreas.  But this has not proved useful for this purpose (at least not without a lot of additional processing).

I then moved on to Linnear Referencing.  Created routes from the stream feature class, and tried to Locate Features Along Routes.  But after leaving it running for hours, I gave up.  Again, even if it worked, I'm not sure I could use the results to achieve what I want.

So my third attempt was to create a Python script which finds all starting nodes of the stream network, then using cursors, traces each route to its terminus by connecting the 'to_node' of one feature to the 'from_node' of the next and builds a single simple line along the way by unioning the geometries of all the connecting line features.  The resulting long line is then inserted (again using cursors) into a 'Routes' feature class.  The idea being that every possible route that water could take is now represented by a single part line feature with no branches.  I'm confident that this could be used to get the results I need - if it ever finished running.  Edit:  This took about 6 days to run, and by then I'd solved it using something similar to the answer provided below.

What is the best way to determine, for each SourceArea, what is the minimum distance to any downstream DangerArea?

Update - new 'Flow Distance' tool:  I've just discovered that there is a new ArcGIS Tool called 'Flow Distance' that looks like it should do exactly what I want.  Unfortunately, it is only available in ArcGIS 10.6 and up, and we are on 10.5.1 - just one version shy of having this tool.


texts: Python GDAL: Write new raster using projection from old

texts:
If I read in a raster image as an array, then make some changes to the values in the array, how do I then save the array as a raster with the same projection information as the original array?

In particular I am conducting processing on some ISIS3 cubes from Mars. These are not projected in any of the nice SetWellKnownGeogCS options. Perhaps this makes my problem somewhat unusual, but I thought it worth documenting my solution all the same.


texts: Where is dataset uuid generated by pyCSW harvest?

texts:
When pyCSW harvests from e.g. a WMS service, it will create a database entry with a uuid.
I wonder, where this uuid comes from.

Is it generated by pyCSW on harvesting or is it inherent to the harvested dataset and only obtained by pyCSW?

If it is generated by pyCSW, how can datasets be updated (based on the same uuid) on subsequent harvesting of the same service?

I didn't find relevant information in the documentation.


texts:
The GRASS plugin toolbar appears incomplete in QGIS 2.14.18.
I can open the mapset but I'm not finding how to load (open) an existing layer.


texts: QGIS-Custom functions not showing in Field Calculator

texts: I'm using Essen and trying to follow Anita Graser's wonderful tutorial Illuminated contours.  I've created the azimuth function and loaded it but there is nothing showing in the expression tab - not even Custom with the default func!  The expression exists as azimuth.py in C:\Users\Max\.qgis2\python\expressions with associated .pyc file.  Restarting QGIS doesn't make any difference


texts:
I am trying to convert a GeoJSON file into an SQLite database and I want to set to INTEGER a specific field in the SQLite final file.

I use the following command line:

  ogr2ogr -f SQLite -mapFieldType my_field=Integer
  final_file.sqlite source_file.geojson

But as a result I get an error:

  ERROR 5: Invalid value for -mapFieldType : my_field=Integer

Without the mapFieldType option it works well.
I used the example shown in the GDAL wiki
https://trac.osgeo.org/gdal/wiki/rfc31_ogr_64


texts: Difference between datum and ellipsoid for geodetic coordinates?

texts:
Looking at the Wikipedia page for the North American Datum, it seems to imply that the datum is completely defined by the ellipsoid, which defines the major/semi-major axes, flattening, etc. For geodetic coordinates, are a datum and an ellipsoid definition the same thing? Can the terms be used interchangeably?

e.g. is NAD83 wholly and completely defined by the ellipsoid GRS80?


texts: Converting coordinates of a bunch of polygons in a shapefile which are in UTM to latitude and longitude

texts:
My exact problem is that I have a file that report vehicular traffic that has 6 columns and two of these columns are latitude and longitude.for example this is the first row and 45. and 9. are latitude and longitude columns.

77c1b49236f442f37531551e2fc25e51    1430352000000000000 45.450824   9.211725    20  0
                         .
                         .
                         .

and on the other hand I have this shapefile of polygons

              id                                           geometry
0           3939_1_1  POLYGON ((504174.2673271392 5003118.268122713,...
1         3939_1_2_1  POLYGON ((508268.9938896392 5003118.268122713,...
2     3939_1_2_2_3_3  POLYGON ((511851.8796318267 5005064.825739901,...
3     3939_1_2_2_3_0  POLYGON ((511340.0388115142 5005064.825739901,...
4       3939_1_2_2_0  POLYGON ((510316.3571708892 5004415.973200838,...

with {'init': 'epsg:32632'}

As you see each polygon has an id I want to check each pair of latitude and longitude and see each one belongs to which polygon and add the column of polygon ids in my first file. So first I read my data file, then I converted the pair of latitude and longitude into a point and I created a geo data frame and then I read the shapefile and tried to check,
This is the piece of code that I have tried

import pandas as pd
import geopandas as gpd
import geopandas.tools
from shapely.geometry import Point, mapping, shape
data_path = "/home/foroogh/PhD/milano-fcd/test"
filename = "2015-04-30-1430352000000000000"
mapDir = "/home/foroogh/PhD/milano-grid"
mapName = "intersection_Milano_W_GRIDIT_NEW.shp"
map_path = os.path.join(mapDir, mapName)
data = os.path.join(data_path, filename)
places = pd.read_csv(data, sep="\t",header = None, names =   ['id','dateandtime', 'latitude', 'longitude', 'vehiclecode', 'velocity'])
places = places[["id", "latitude", "longitude"]]
places["geometry"] = places.apply(lambda row: Point(row["longitude"],  row["latitude"]), axis = 1)
del(places["latitude"], places["longitude"])
places = gpd.GeoDataFrame(places, geometry = "geometry")
places.crs = {'init': 'epsg:32632'}
milano = gpd.GeoDataFrame.from_file(map_path)
print(milano)
milano = milano[["id", "geometry"]]
print(milano)
result = gpd.tools.sjoin(places, milano, how = "left")
print(result.head())

but I do not get the desired result and the problem is that my shapefile is in utm and my points in geo dataframe are in lat and lon and I did not know which one should be converted and how?


texts: Counting points and adding to column with shapely

texts:
I have a small dataframe and created new column. I would like to add to every record in column 'vertex_counter' value of n. I have the problem because it is doesn't work and I have in every record in column 'vertex_counter' value 0.
What is the problem?
df_budynki3 = df_budynki
df_budynki3['vertex_counter'] = 0

for index_budynku, bud_row in df_budynki3.iterrows():
    coo = shapely.wkt.loads(bud_row['classification.geometry']) <-- create polygon from string
    n = len(coo.exterior.coords) <-- count number of points in polygon
    bud_row['vertex_counter'] += n


texts: Geoserver JDBC User Group Services Problem

texts:
I use the Geoserver 2.11.4 version in Windows.I am adding the JDBC service to the PostgreSQL database using JDBC Users and Group Services. Next, I add a Provider for JDBC that I created using the JDBC Authentication Provider job queue. After doing this, the tables are created in the Database and I can connect as admin with postgres user. When I add user from JDBC service, user information is written in the database.
PostgreSQL view

but I can not login to Geoserver with the user1 user I added. I can login with the user I created when I added the user from his user adding field even though it can not be in the table in Postresql. But I can not assign a rule to this user. How can i this problem.


texts: Black Screen problem during boot of osgeo-live 10.5 with VirtualBox on Windows10

texts:
Has anybody been able to use the osgeo-live virtual machine disk with VirtualBox?

I follow the instructions from the osgeo website but the machine will not boot.

I am stuck at the GNU GRUB selection. All boot options are a black screen with an empty cursor that does not take any input.

My default params are:

The advanced options for Ubuntu boots with the message 'Initialising Ram disk' before going dark on me while the default options starts black.


texts:
I am working on a paper related to the land cover changes between 1984/85 and 2010/11 in my country and I am searching for Landsat 5 images from December 1984(2010) to March 1985(2011). Then I average band-by-band using raster calculator to build a single image for each period to evaluate land cover changes. Before averaging, each band will be corrected (ToA correction) and the image averaged will be corrected geographically.

Is this procedure correct?


texts: Limit ST_DWithin() to first candidate found

texts:
Im trying to create a table with distances calculated from points to lines.
For that im using the following SQL adapted from BOSTON gis website.

I understood the function seems to search for all candidates within the specified distance, sort them and pick the closest one

Happens that i would like to limit the search when it founds the 1st match within   the specified distance. (hoping for faster results)

Should this be done diferently ?

DROP table if exists filter_distance.table_distances_@country;
CREATE table filter_distance.table_distances_@country as (
SELECT DISTINCT ON(g1.gid)  g1.gid As gid, g2.gid As id_lines , round(cast(ST_DistanceSphere(g1.geom,g2.geom) as Numeric), 2) as distance
    FROM filter_distance.points_@country As g1, filter_distance.polys_lines_@country As g2
    WHERE g1.gid <> g2.gid AND ST_DWithin(g1.geom, g2.geom, 0.00500)
    ORDER BY g1.gid, ST_Distance(g1.geom,g2.geom) asc);


texts: Importing 3D Analyst functions with ArcPy?

texts:
This really seems like I should be able to do this, but I am running into some problems. I want to import 3D Analyst with python, in the same way I import Spatial Analyst, but it returns an error.

For Spatial Analyst I use the following:

from arcpy.sa import *

This allows me to use the spatial analyst functions in python scripts outside of ArcMap. I want to do the same with 3D analyst, but I am running into an error.

I really feel like it should just be:

from arcpy.3D import *

I also tried arcpy.3d, but either way I get the error:

Parsing error <type 'exceptions.SyntaxError'>: invalid syntax (line 1)


texts: ERROR 000539 from ModelBuilder with Raster Calculator?

texts:
I'm trying to do a math operation in raster calculator using the rasters from other two raster calculators.

My calculation is a simple NDVI formula:

(rastcalc4 - rastcalc3)/(rastcalc4 + rastcalc3)

It works well in most of the cases, but, sometimes it fails and this error shows up:

  ERROR 000539: Error running expression: rcexec()

  Traceback (most recent call last):

  File "", line 1, in

  File "", line 8, in rcexec

  RuntimeError: ERROR 999998: Unexpected Error.

  Failed to execute (Raster Calculator (3)).

I don't understand what that means, how can i fix it?


texts: Understanding difference between area calculation when using zonal statistics tool and identify button?

texts:
I am working in ArcGIS Desktop 10.5 and I am a beginner.

I have very basic, non-overlapping polygons over a raster that represents bathymetry data.  I need to figure out various statistics about each polygon (like range) so I figured I would use the Zonal Statistics tool.  I also needed the area of the features. So I did this and things were fine.

But then I happened to notice that the areas I got for my polygons were different from the "Shape_area" calculations that I see when I use the "identify" tool (the one in a blue box with an 'i' on it) and click on the polygon.

Which area is correct?


texts: Does the strength of GPS weaken in rural areas? For example, based on how the satellites work, would the strength in the centre of Australia differ from a populated city in Europe?


texts: In QGIS 2.16.1 and 2.16.2 no Raster Calculator available anymore

texts:
In all my older versions the raster calculator is easy to find in the raster menu - however in both 2.16. versions, it is not visible. As I understand this does not have anything to do if GdalTools is active or not.

I heard from one other user that she can see the Raster calculator in 2.12.2.

My setup is on a windows 10 PC.

Anyone else with the same problem?


texts: Turning off only some of layers specified using ArcPy?

texts:
I only know how to turn off the entire current layer. But I want off only some of the layers specified

#turn all layers in current map doc on:
print "Button2--turn on layers"
mxd = arcpy.mapping.MapDocument("Current")
list_layers = arcpy.mapping.ListDataFrames()
for layer in list_layers:
    layer.visible = True


texts:
I have a shapefile displaying the borders of several city units. One of its polygons, however, had free space consisting of four vertices in it. I created another polygon inside the free space and then used Dissolve to close the free space. The result was the following layer which still contains some errors (see the close-ups below).

Is the best way to resolve this in QGIS 2.2 using some kind of v.clean?

The entire data seems to be corrupted. Using Dissolve on the entire shape results in what the fourth screenshot shows.


texts:
I have two raster layers with different resolutions that I want to join.

One has higher resolution (transparent yellow) and the other raster has lower resolution but bigger extent (whole earth) and has information about different classes (drawn in different colors here). The resulting raster should have the higher resolution and extent of the raster drawn in yellow here but should be joined with the other raster, e.g. containing the information of what class it was laying within.


texts: Why does ModelBuilder with kriging give ERROR 001000 Z value field does not exist?

texts:
I've built a model which interpolates series of data by IDW method and it works well. But when I tried to change method of interpolation to kriging I encountered some problems with using Z value field as a variable.
Here are some pictures that may help to explain my problem:
how the data look like: http://imageshack.us/a/img22/4167/l4th.jpg IDW version: http://imageshack.us/a/img221/1879/qegc.png
I'm using 'For' iterator and variable for Z value field 'a%Value%' so every iteration is using next data field as Z values. And values of the iterator are preconditions to interpolation. It's is basically the same model in both versions but I don't know why the variable isn't working with kriging. When I open the Kriging tool I've got this error in Z value field:

  ERROR 001000 Z value field: Field a%Value% does not exist

but in the IDW version it works well. Any ideas?


texts: How open Web AppBuilder widget programmatically?

texts:
I'm trying open a widget from other widget:

_open: function () {

    var config = WidgetManager.getInstance().appConfig.getConfigElementById("widgets_IdentifierRed_Widget");
    WidgetManager.getInstance().loadWidget(config).then(function(widget){
      WidgetManager.getInstance().openWidget(widget.id);
      console.log(widget);
    });            }

First I load it, then open it. But it doesn't work.

How can I open a widget from other widget?


texts:
I try to write a PyQGIS python script to work standalone in the python IDE without GUI, out of QGIS interface following the oficial tutorial.

The tutorial for pyqgis in standalone scripts says:

    from qgis.core import *

    # supply path to qgis install location
    QgsApplication.setPrefixPath("/usr", True)#for linux

    # create a reference to the QgsApplication, setting the
    # second argument to False disables the GUI
    qgs = QgsApplication([], False)

    # load providers
    qgs.initQgis()

    # Write your code here to load some layers, use processing algorithms, etc.

    # When your script is complete, call exitQgis() to remove the provider and
    # layer registries from memory
    qgs.exitQgis()

If I write print 'hello world' it shows that message succesfully

from qgis.core import *
# supply path to qgis install location
QgsApplication.setPrefixPath("/usr", True)#for linux

# create a reference to the QgsApplication, setting the
# second argument to False disables the GUI
qgs = QgsApplication([], False)
print 'hello world'
# load providers
qgs.initQgis()

# Write your code here to load some layers, use processing algorithms, etc.

# When your script is complete, call exitQgis() to remove the provider and
# layer registries from memory
qgs.exitQgis()

The tutorial says # Write your code here to load some layers, use processing algorithms, etc. after the line qgs.initQgis() and before qgs.exitQgis().

But if i add lines for processing plugins like :

qgs.initQgis()
sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
import processing
Processing.initialize()
processing.alghelp("saga:slopeaspectcurvature")
qgs.exitQgis()

These commands don't work, they don't show me anything.

If i write that command lines in other places in the script it dont show me anything, not even my hello world.

If i write sys.path in the beginning of the script it doesn't show me anything again.

Do you have any idea how to write the right code to get access to the processing tools with standalone pyhon script?

On the python console inside the QGIS import processing and processing.alghelp("saga:slopeaspectcurvature") work fine.


texts: Finding all records related to selected records in ArcGIS Desktop using ModelBuilder?

texts:
I have two layers that are related. My current manual process is as follow;

a.  Add the layer 1 to the table of contents

b.  Open the attribute table then Switch selection to select all layer 1 records

c.  Click relates and select layer 2
This will select all layer 2 records related to the layer 1

d.  Switch selection
This will display all layer 2 records not related to the layer 1

e.  Export the table to excel

I want to automate this using ModelBuilder.

This is where I am using ModelBuilder;

Layer 1 > Select Layer By Attribute (Switch selection) > Calculate Field (I require a python script that is able select all related layer 2 records) > Select Layer By Attribute (Switch Selection) > Export to Excel


texts: Google Engine - NDVI time series by region - LineChart not connecting dots

texts:
I created a stack of NDVI images, and also have several polygon. When I plot the time series, points show up but only a few connects dates to each other.

I am using the ui.Chart.image.seriesByRegion()

https://code.earthengine.google.com/?scriptPath=users%2Fdelaralfonso%2FAPPs%3AExample

Any thoughts how to fix this?


texts: Assigning coordinate system to ungeoreferenced CAD Drawing while maintaining CAD shape and size using ArcMap

texts:
I am using ArcMap 10.8.1 for this project. I am attempting to bring a small town's CAD map of their streets and buildings into ArcMap and then will work on converting some of this data to shapefiles.
The CAD Drawing doesn't seem to have any sort of coordinate system, but I can use the measure tool within CAD and everything seems accurate there (the scale in the legend matches up with the measuring tool).
However, as soon as I bring it into ArcMap the drawing gets massive and is over 10x bigger than it should be. The scale bar should be 2000ft but when I measure it in ArcMap it is close to 30,000!
I am using the correct projected coordinate system from what I can tell (NAD1983 Illinois East) and the GCS is also NAD 1983, so I am unsure of what else I can do to match it up. I've tried playing with these a bit but with no luck.
I've never done something like this before, so I'm 100% sure I'm just missing a step here. I was able to move the drawing on top of a reference point I made using the BaseMap, but that's not much help if the scale is completely off.
Is there anything I am missing here?


texts: Retaining co-ordinate values with 6 significant figures when converting shapefile into Excel spreadsheet?

texts:
Using ArcGIS-desktop, I converted a shapefile, into an Excel table through the conversion tool Table to Excel, however my co-ordinates previously were specific to the 1x10-6 (i.e. 9.078627) but in the new Excel spreadsheet they have been converted into whole numbers (i.e. 9).

Is there a method to retain the original specificity of the data points when converting?


texts: Select lines on top of lines in same shapefile with ArcGIS

texts:
I have a street centerline shape file where segments for bi-directional roadways are drawn as two stacked lines. The top segment for one direction of travel, the bottom segment for the opposite direction of travel.
I need to separate these stacked lines into two different shape files.  Anyone know how to select just the top set of lines so I can export them to a different file? I can just use switch selection to get the second set of lines then.
Unfortunately none of the table attributes represent the top versus bottom segment, or direction of travel.  The geometry of the top and bottom segments are the same.
I looked at commands such as Find Identical, but that still selects both the top and bottom at the same time.
Looking for a solution in ArcGIS Desktop 10.6.1 Advanced or ArcGIS Pro 2.6.2


texts:
I am attempting to calculate albedo using NAIP imagery 4 bands (visible-NIR).

Zhao et al 2015 used quick bird imagery to calculate albedo using the following equation:

Albedo = -1.694(band1(blue))+2.749(band2(Green))+0.288(NIR).

Can I use these same coefficients that Zhao et al used?

Quickbird has band designations as:

Blue: 430 - 545 nm
Green: 466 - 620 nm
Red: 590 - 710 nm
Near-IR: 715 - 918 nm

and NAIP has band designations as:

Blue: 420-492 nm
Green 533-587 nm
Red 604-664 nm
Nir: 833-920 nm


texts: Identify spatially contiguous clusters in raster data using kmeans

texts:
I would like to cluster the cells of a raster object into k contiguous regions using kmeans. The number of regions, k, is known. Each cell has various geographical attributes, such as temperature, precipitation, elevation etc. And I need R to sort pixels into k groups (regions) using these cell values. With stats::kmeans() that is a fairly simple exercise. Unfortunately, this method does not create spatially contiguous clusters. Instead, each group consists of pixels spread all over the grid.

I expected this to be a common problem, but could not find any R function that solves it. There should be some package out there that can perform clustering under the restriction of spatial contiguity.

The method of choice here is kmeans because I know how many groups I need. Any other technique from the field of unsupervised learning that allows me to a priori set this quantity, is of course also welcome.

I asked this question on Stack Overflow before and received a recommendation to also post it here, as GIS might be a better fit than SO.


texts: Cannot open zip, compressed, file on external drive with GeoPandas or Fiona, it opens on C drive

texts:
I am working with Windows 10 and GeoPandas version '0.4.1' which utilizes Fiona version '1.8.11'

I am trying to open a shapefile from a compressed folder (as described here https://geopandas.org/io.html) on an external hard drive.

my code:

gpd.read_file(r"zip:///ecan_data_delete/Waimakariri/Waimakariri/Groundwater/Numerical GW model/input_n_soil_conc.zip!input_n_conc_2018-11-24.shp")

returns the following error only when the file is located on an external harddrive (e.g. zip file path "D:\ecan_data_delete\Waimakariri\Waimakariri\Groundwater\Numerical GW model\input_n_soil_conc.zip"):

fiona.errors.DriverError: '/vsizip//ecan_data_delete/Waimakariri/Waimakariri/Groundwater/Numerical GW model/input_n_soil_conc.zip/input_n_conc_2018-11-24.shp' does not exist in the file system, and is not recognized as a supported dataset name.

when I make a copy of the zip file on my local drive (e.g. create "C:\ecan_data_delete\Waimakariri\Waimakariri\Groundwater\Numerical GW model\input_n_soil_conc.zip")
it runs perfectly, which means that fiona is assuming that the file (and therefore I believe the virtual file) is on the C drive.

Is it possible to fix this?  I have spent some time looking through the documentation, but it is pretty light here (and/or above my understanding).  I know that I could simply move the zipped file over to my C drive, but I am running through a bunch of shape files zipped up with some big modelling files on a hard drive and making some quick figures to see what is there, so moving it is not a very efficient solution.


texts: Why is my attribute table edit button greyed out?

texts: I'm very new to GIS technology and have been trying to learn qGIS through downloaded tutorials. My qGIS is the 2.14.21 version. It was going well until the instructions told me to edit the attribute table of one of my layers (a sqlite layer with a polygon geometry type).  The edit button is 'greyed' out, not allowing me to access this function. Why is this? I have tried various ways of selecting the data, restarting the program, and accessing the table through different buttons. The 'add a new field' button is also greyed, though the field calculator is functional. My online research hasn't come up with answers yet.


texts: Styling KML through LIBKML layer creation options

texts:
I am exporting a layer from PostgreSQL to KML using ogr2ogr.

I get a workable KML with a series of features (Placemarks). Here is a simplified example of a feature. (I took out the long list of polygon points to simplify the length of the geometry).

<Placemark>
    <name>FirstFeature</name>
    <description>A map feature</description>
    <Style>
      <LineStyle>
          <color>ff0000ff</color>
      </LineStyle>
      <PolyStyle>
          <fill>0</fill>
      </PolyStyle>
    </Style>
    <ExtendedData>
        <SimpleData name="MapUnit">myunit</SimpleData>
        <SimpleData name="hexcolor">ffe9e9e9</SimpleData>
    </ExtendedData>
    <MultiGeometry>
         <Polygon><outerBoundaryIs><LinearRing>
             <coordinates>[array of coordinates]</coordinates>
         </LinearRing></outerBoundaryIs></Polygon>
    </MultiGeometry>
</Placemark>

I would like to set the layer creation options such that it will refer to the 8bit hex-code provided by the hexcolor attribute. By default, ogr2ogr seems to be styling all polygon features as red lines with transparent fill (see above <Style></Style>). I instead want the style such that the <Style></Style> section becomes...

<Style>
    <LineStyle>
        <color>ff0000ff</color>
    </LineStyle>
    <PolyStyle>
        <color>ffe9e9e9</color> <!--Taken from the hexcolor attribute-->
        <fill>1</fill>
    </PolyStyle>
</Style>


texts:
How do I display shape files with a google-like interface in a web application?

I found a tutorial at http://spatialhorizons.com/2007/05/19/mapserver-1-10-minute-tutorial/
which uses an earlier version of map server (the tutorial uses v2.2.6 and I am using v3.0.6) and followed it with no problems up to the second to last point.  I know I have set up map server correctly since when I navigate to localhost in my browser I see the welcome screen.  However, when I select the ka-map sample interface link I get the error: "The requested URL /kamap was not found on this server."

Could someone please advise about where I am going wrong?


texts: How to make infobox contents appear on mouseclick instead of mouseover

texts:
I'm using CartoJS to put one of my CartoDB maps on a div in an HTML file. When a user mouses over a dot on the map, info appears in the div next to the map container div.

I want to change this so the info appears when a user clicks (or taps) a point.

Here is a snippet of my JS code for making the CARTO layer.

            cartodb.createLayer(map, layerURL, {infowindow: false, tooltip: false})
                .addTo(map)
                .done(function(layer){
                    var infoBox=    layer.leafletMap.viz.addOverlay({
                        type: "infobox",
                        layer: layer,
                        template: infoBoxTemplate,
                        position: "top|right"
                    });

                    // Add infoBox to map object
                    $("#infobox-wrapper").append(infoBox.render().el);
                });

How do I change this so info appears in #infobox-wrapper on click instead of on hover?


texts: How to calculate the area of polygon features within another polygon feature layer

texts:
As shown below, I would like to calculate the percentage of each red feature (landscape zones) taken up by a green feature (forest zone). So I need the summed percentage of forest zone for each landscape zone feature. I have a unique ID field each of the landscape zone features. I think I would need to dissolve all the forest zones that fall within each red feature into single features. However, I am unsure how to do this - would it be a spatial join? I believe I could then simply create a new field and divide the forest zone areas by the landscape zone areas to get a percentage of forest zone area for each landscape zone feature. In total I have over 800 landscape zones to do this for.

I have both ArcGIS and QGIS and am looking for the easiest way to do this.


texts: Field calculator Python else statement question

texts:
The codeblock below is run on Field1 and calculates Field1 to equal Field2 if they are not equal.  Else it returns Field1.  My question is, does the else statement actually recalculate Field1 to equal itself or does it leave Field1 as is.  Basically if Field1 is equal to Field2, I don't want Field1 to recalculate to itself.

def update(Field1,Field2):
  if Field1 != Field2:
    return Field2
  else:
    return Field1


texts: QGIS Unknown CRS from geotiff created with python GDAL 2.4

texts:
I've been creating geotiffs using python gdal, here is a minimal example of my code:

from osgeo import gdal, osr

driver = gdal.GetDriverByName('GTiff')

dataset = driver.Create("test.tif", 20, 10, 1, gdal.GDT_Float32)
dataset.SetGeoTransform((500000, 5, 0, 100000, 0, -5))

srs = osr.SpatialReference()
srs.ImportFromEPSG(27700)
dataset.SetProjection(srs.ExportToWkt())

dataset = None

I can then attempt to view this geotiff in QGIS (3.12.1).  If I use the above code using GDAL 3.0, then the CRS is recognised, that is, in the layer properties I see:

EPSG:27700 - OSGB 1936 / British National Grid - Projected

as I expected.  However, I had been using GDAL 2.4, and using the above code gives:

Unknown CRS: BOUNDCRS[SOURCECRS[PROJCRS["OSGB 1936 / British National Grid",BASEGEOGCRS["OSGB 1936",DATUM["OSGB 1936",ELLIPSOID["Airy 1830",6377563.396,299.3249646,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4277]],CONVERSION["British National Grid",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",49,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-2,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996012717,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",400000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",-100000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["unknown"],AREA["UK - Britain and UKCS 49�46'N to 61�01'N, 7�33'W to 3�33'E"],BBOX[49.75,-9.2,61.14,2.88]],ID["EPSG",27700]]],TARGETCRS[GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["unknown"],AREA["World"],BBOX[-90,-180,90,180]],ID["EPSG",4326]]],ABRIDGEDTRANSFORMATION["Transformation to WGS84",METHOD["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",446.448,ID["EPSG",8605]],PARAMETER["Y-axis translation",-125.157,ID["EPSG",8606]],PARAMETER["Z-axis translation",542.06,ID["EPSG",8607]],PARAMETER["X-axis rotation",0.15,ID["EPSG",8608]],PARAMETER["Y-axis rotation",0.247,ID["EPSG",8609]],PARAMETER["Z-axis rotation",0.842,ID["EPSG",8610]],PARAMETER["Scale difference",0.999979511,ID["EPSG",8611]]]] - Projected

I would appreciate if someone could point me to the reason why there is a difference.  And also, why QGIS does not recognise the version from GDAL 2.4. Thanks in advance.

Just to add to the above. I tried to see if there was a difference in the definition of the spatial refs in the geotiffs. I used dataset.GetSpatialRef().ExportToProj4() on the datasets above and they yielded:

'+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs' # gdal 3
'+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs' # gdal 2.4

Note that the only difference is the presence of towgs84 in the gdal 2.4 version. Is this relevant to the problem I'm seeing?


texts: Can Spatial Joins be used to perform One-to-Many Closest joins?

texts:
My situation appears rather simple, but it has given me headaches more than once already. I have a set of points and a set of polylines. I would like to use a spatial join to join the attributes of the closest line to the point data. However, there are many lines that have spatially coincident counterparts. That is to say, duplicate geometry exists that allows more than one line to occupy the exact (emphasis) space as another line. However, when performing a One-to-Many join for Closest feature in a Spatial Join, only one line is chosen. This does not seem to be proper behavior, but perhaps I am mistaking the terminology? If so, is there a better method to create an output table that creates duplicate point features with the attributes of all lines that are "closest?"
Pictures for clarity:

As you can see, there are 5 line features that are coincident at this location, but a one-to-many spatial join for closest feature only picks one. Coincidentally (that word never gets old), it is the first feature identified when the link is clicked with the Identify tool. They appear to list in order of ObjectID, but this may just be a... well, you know.
Any ideas?


texts: Cant import and show a vector in PyQGIS?

texts:
This is what I have done.

from qgis.core import *
layer = iface.addVectorLayer(r"C:\\Users\\username\\Desktop\\coastline.shp","coast","ogr")

or

from qgis.core import *
layer = QgsVectorLayer(r"C:\\Users\\username\\Desktop\\coastline.shp","coast","ogr")

Basically it doesn't even appear in the Layers section of QGIS.

The first one says layer is not valid
and the second it doesn't cause any problems but it isn't in the layers at all of Qgis?


texts: PyQGIS - categorized renderer classification - create category with 'no value' or 'any other value'

texts:
The code snipped mentioned in Automatic pyqgis categorized renderer classification and Apply symbol to each feature (Categorized symbol) automatically create categorized renderers and determine categories using PyQGIS.

I (using QGIS 2.18.18 LTR) have a layer with uniqueValues for a field pre_cluster_id like [0, 1, 2, 4, 5, 6, 7, 8, 9, NULL] for which the above referred to approaches correctly create a category with value 0 and Legend 0, but incorrectly (from my understanding) create a category with value 0 and legend NULL hiding the first one, screenshot:

Created with the Classify Button, there is a category with no value, covering all except [0, ..., 9] (including NULL):

Hence, question: Using PyQGIS, how can I create a QgsRendererCategoryV2() with no value, i.e. all other values (or a category for NULL value, respectively)?


texts:
I am trying to show dojo GFXLayer of dojox.geo.openklayers on my map layers , the problem is that my map is not showing townLayer but map is showing default OSM Layer and GFX layer.What I am Doing Wrong?
My Code Is

var dojoMap = new Map("map");
var map = dojoMap.getOLMap();
var townLayer = new OpenLayers.Layer.WMS("Towns","http://"localhost":6666/geoserver/wms?",
   {layers: 'Layer:towns', format: 'image/png', transparent: true, info_format: 'text/html',
    gtype:'polygon'});

map.addLayer(townLayer);
var layer = new GfxLayer();
var p = new Point({x:ny.longitude, y:ny.latitude});
var f = new GeometryFeature(p);
f.setFill([ 0, 128, 128 ]);
f.setStroke([ 0, 0, 0 ]);
f.setShapeProperties({r : 20});
layer.addFeature(f);
dojoMap.addLayer(layer);


texts: GEE - Get date of ImageCollection mosaic

texts:
I have created mosaics by date from an ImageCollection and now want to label each mosaic with its corresponding date. This is my function to label the ImageCollection:

var date_label_pos = ee.Geometry.Point([32.55, -15.42]);
var scale = 200;

function annotate(collection) {
  return collection.map(function(img) {
    // define text properties
    var props = {
      fontFamily: 'Calibri',
      fontStyle: 'bold',
      textColor: '1D4180',
      outlineColor: 'ffffff',
      outlineWidth: 2.5,
      outlineOpacity: 0.6
    };
    var center = ee.List(date_label_pos.coordinates());
    var pos = translate(center, 0, 0);
    var s = ee.String(' DATE SHOULD BE HERE!');
    var date_label = Text.draw(s, pos, scale * 1.5, props);
    // merge results
    var image = img.visualize({forceRgbOutput: true, min: 500, max: 1800, bands: ['B4', 'B3', 'B2']});
    return ee.ImageCollection([image, date_label]).mosaic();
  });
}

The full code is HERE. What I'm missing is the part where I can extract the date from the image mosaic. img.date() does not work.

To provide a little bit more info: I needed to set a new property date to the mosaics, which I can see when I print the ImageCollection to the console. I can also access that date outside of the annotate function using:

print(ee.Image(mosaics.toList(mosaics.size()).get(0)).get('date'));

where mosaics is my ImageCollection. But within the function, it seems that img is something different and I can't figure out what it is and how I can access the date property I set.


texts: ERROR: more than one row returned by a subquery used as an expressionSQL state: 21000 Why and how it happens?

texts:
I want to update an empty column with a join :
UPDATE icpe_modifie
SET region =
(SELECT nom_reg_m FROM icpe_modifie as icpe
JOIN regions as reg
ON ST_WITHIN(icpe.geom, reg.geom))

icpe_modifie is a point table and regions is a polygonal table.
I want to add the name of the region where is located the point.
But it returns :
ERROR:  more than one row returned by a subquery used as an expressionSQL state: 21000

So I tried :
SELECT nom_reg_m FROM icpe_modifie as icpe
JOIN regions as reg
ON ST_WITHIN(icpe.geom, reg.geom)

This second query is well executed and no error is raised.
The current solution I chose is to LIMIT 1 :
UPDATE icpe_modifie
SET region =
(SELECT nom_reg_m FROM icpe_modifie as icpe
JOIN regions as reg
ON ST_WITHIN(icpe.geom, reg.geom) LIMIT 1)

But I would like to better understand the error statement, to it means that join returns for at least one row several items.
I do not understand how a point could be in two different areas. And I also would like to check the points for which this is happening and have no idea to find these trouble points.


texts: Error [37] : Unable to calculate. Metadata error - SCP plugin in QGIS - Preprocessing landsat 8

texts:
I am trying to preprocess Landsat 8 data using the SCP plugin - all the bands are processed except for bands 1 and 9. The metadata is read correctly in the table but on clicking "Run" I hit upon the error.

Error [37] : Unable to calculate. Metadata error - SCP plugin in QGIS
Preprocessing Landsat 8

I'm using QGIS 3.16.2.


texts: Does ArcGIS Online include High Speed Rail in its network analysis?

texts:
I'd like to run an OD Matrix using network analyst in Europe and Asia to study high speed rail.

Does ArcGIS online include HSR (such as the AVE in Spain or the TGV in France) in its analysis?


texts: Convert 3D-Buildings as MultiPolygonZ to 2D shapefile with height feature

texts:
As of March 2021 Switzerland provides the swiss 3D building geodataset for free, it is available as dxf or gdb, the latter contains MultiPolygonZ geometries (buildings also with complex roof geometries). I want to convert these to a 2D shapefile, where building height is a field. I see that the derived geometry in the attribute table contains an interpolated Z and a Next support point Z (in german: Nächster Stützpunkt Z).
Is there a tool or process in QGIS to reduce the 3D dataset to 2D? Or would it be easier to iterate over every MultiPolygonZ in geopandas, keep the extents and average the uppermost z support point values?


texts: QGIS Processing algorithms have encoding problems

texts:
I use QGIS 2.14.1.

I add vector layer from UTF-8 encoded geojson, which is read fine and attribute table show correct characters.

If I use QGIS algorithm on this layer (i.e. clip, voronoi, etc) I get new layer which attribute table is full of question marks instead correct characters.

If I use SAGA algorithm (i.e. clip) characters in newly created layer are fine.

Is there some setting I'm missing?


texts:
I'm sure this shouldn't be as tricky as I'm finding it but does anyone know how you can find the 'remoteness' of a point compared to multiple other points in QGIS?

I have a shapefile containing 1000's of points, and I want to find out which one of them is the most remote. What I mean by remote is:

"the state of being distant from something else"


texts: Insights in AGOL - "We couldn't display all features" error

texts:
I have a file geodatabase with a single feature class (polygon) which I have zipped and uploaded to my organisation's AGOL account and created a hosted feature layer.

There are just 11 polygons in my feature class.

I can view the data in a web map.

I can view the data in ArcMap.

I have tried adding it to Insights but it won't display all polygons when I add it to a map in Insights.  Some polygons will only display if I zoom in really close.  Some polygons won't display at all.  The 11 polygons encompass all the water areas around the UK, so they aren't small polygons.

If I open the dataset as a table in Insights, I can see all the records.

The same dataset was working a few weeks ago when I last used it in Insights.  The dataset has not been edited since.

The error I get:

If I click for more information, it just tells me how many features are being displayed.

The data is in wgs84.

Things I have tried.

Chrome (v66) and Internet Explorer 11.
Cleared cookies & cache.
In ArcMap, checked for geometry errors - there are none.
Saved the feature class to a new fgdb and added it to AGOL.  Same error.
Exported the feature class to a shp file, zipped this and added it to AGOL.  Same error.
Deleted the .prj file for the shp file and run the "Define Projection" tool in ArcMap and set it to wgs84.  Same error.
Re-projected the dataset to Web Mercator. Same error.
Checked the extent of the dataset.  This is correct.  As shown by the automatic thumbnail created in AGOL. I have double checked this by zooming to extent in ArcMap and AGOL.

Re-created the dataset from scratch.  Same error.
In the settings tab for the hosted feature layer in AGOL, I've optimised the layer drawing.  Same error.
Published as a service straight from ArcMap.  Same error.
Generalised the polygons in ArcMap and then repaired geometry to fix any self intersections caused by this process.  Have tried a 10m, 30m and a 100m generalisation.  At 100m, the vertices on the most complex polygon reduced from 50,000 to 19,000. However, there is no correlation between the polygons with the highest number of vertices and the polygons not displaying. i.e. The polygons with the fewest vertices (1233 pre generalisation, 469 post generalisation) is one of the ones not displaying.

I am at a total loss for what to do next.  Any ideas?

Is this just a glitch with Insights?  It's just confusing as the same dataset was working a few weeks ago.


texts: What's going wrong when I reproject raster in EPSG:32115 in QGIS?

texts: I wish to make a QGIS project with an orthoimagery layer and a vector layer of GPS points. I'm using QGIS 2.8.2, 64 bit, in Windows 7 on a Lenova t430 laptop. Raster is in EPSG: 32115 (mercator, NAD83, Eastern New York): I add that layer first. When I add a *.gpx layer, which is in EPSG4326 (Lat/Long) WGS84, the points display 1000 km south from where they should be. When I add the layers in the opposite order, the raster layer ends up way up in Canada somewhere. No amount of reprojecting of either layer seems to be helping. Do you have any suggestions?


texts: Defining a WHERE clause from a subquery for a query to create a handy view

texts:
It is about my implementation from an answer of this question.
It is more an postgresql issue than a GIS issue, whereas it is about geometries and spatial queries.

I want to create a view from the query below, whereby I can select the result by an ID (that defines a linestring).
My problem is, that I had no success with coding the query in that way, that for every linestring only one newpoint is created that is half the length away from the linestring (see last subquery 'newdistance').
I tried to SELECT DISTINCT ON (id)... but that failed.

The linestring table:

id smallint NOT NULL,
geom geometry(linestring, 31468)

I have the same problem in another query, to "extract" a WHERE clause (here in the first part of the query) from an subquery.

My goal is to have a view newpoint (from the query below) whereby I can query like this

SELECT * FROM newpoint WHERE id=1

The query:

WITH line AS
    (SELECT
        geom
    FROM linestring
        WHERE id=1),
points AS
    (SELECT
        ST_StartPoint(geom) AS A,
        ST_EndPoint(geom) AS B
    FROM line),
azimuth AS
    (SELECT
        ST_Azimuth(A,B) AS azimuth
    FROM points),
distance AS
    (SELECT
        ST_Distance(A,B) AS length
    FROM points),
newdistance AS
    (SELECT
        length + (length * 1/2) AS newlength
    FROM distance)

SELECT
    ST_SetSRID(ST_Translate(A, sin(azimuth) * newlength, cos(azimuth) * newlength), 31468) AS newpoint
FROM newdistance, azimuth, points

Setting the WHERE clause id=1 at the end of the query didn't succeed.

EDIT

When I use

SELECT DISTINCT ON (geom) ...

in the last part of the query it works ONLY with 2 objects (linestrings). With more than 2 linestrings it doesn't work.

Here is a picture, how does it has to look like.


texts:
I try to add a buffer on a set of lines via the buffer attribute of geopandas.
But by testing on a simple case (see code below) I realize that the size of the buffer is not kept according to the points latitude.
The buffer build at point A is well 2° but at point B is only 0.35° (checked via QGIS).
Could someone please explain to me how this buffer works under geopandas?
And if you have a solution to generate a buffer that is robust to the change of latitude I'm interested!
import geopandas as gpd

from shapely.geometry import LineString, Point

# 2 points creation
A = Point(0.0, 0.0)
B = Point(0.0, 80.0)

# Linestring creation
lines_geom = LineString([A, B])
crs = {'init': 'epsg:4326'} # WGS84 CRS
lines = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[lines_geom])

# Buffer generation
lines_buffer_1deg = lines.buffer(1.0, resolution=1)
print(lines_buffer_1deg)

lines_buffer_1deg.to_file("test.shp")

As proposed in the comment I also tried with a projected CRS, EPSG:3395 world mercator with the same result (see code below and caption from QGIS).
import geopandas as gpd

from shapely.geometry import mapping, Polygon, LineString, Point

# 2 points creation
A = Point(0.0,0.0)
B = Point(0.0,9000000.0)

# Linestring creation
lines_geom = LineString([A,B])
crs = {'init': 'epsg:3395'}
lines = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[lines_geom])

# Buffer generation
lines_buffer_100km = lines.buffer(100000.0,resolution=1)
print(lines_buffer_100km)

lines_buffer_100km.to_file("test.shp")

EPSG:3395 around 0° latitude:

EPSG:3395 around 62° latitude:


texts: Want data from multiple rasters, but only for given points

texts: Working in ArcMap 10, I have a number of raster layers that are derived from a DEM (terrain attributes, some zonal statistics layers, flow accumulations).  I have a point layer that contains about 25 points that correspond to features I am interested in on the landscape.  As an early step, I would like to pull data from each of the raster layers, but I am only interested in the raster values that correspond to each of the points from my point shapefile.  Essentially I want to make a big table that lists many attributes (elevations, stream power index, flow accumulation, etc.) for each of the 25 points so that I can try to determine if there are relevant trends amongst the 25 points.  How can I do this?


texts: Finding relationship/correlation between two rasters?

texts:
I have two rasters: one with categorical values (vegetation types) and one with a continuous variable (lets say rainfall in mm), and I want to find the relationship.

From the start I thought it would be easy to do it simply in ArcMap or just export the data to Rstudio (or similar) to do statistical analysis for significance value etc.

I have not found a way yet, I  was only able to get the mean rainfall for each category using zonal statistics.

Any ideas?


texts: How to create spatial tag clouds or "tag maps"?

texts:
I'm looking for implementations of spatial tag clouds or "tag maps". Some papers I found talk about a Yahoo Tag Maps API, but I couldn't find any such thing.

I'd prefer OpenLayers-based solutions, but others would be interesting too.

Edit: Example images can be found at http://www.soi.city.ac.uk/~sbbb717/go2/iv_figs/, e.g:


texts:
I'm creating my layer like this:

crs = core.QgsCoordinateReferenceSystem(4326, core.QgsCoordinateReferenceSystem.EpsgCrsId)
uri = "Polygon"
uri += '?crs=' + crs.authid()
uri += '&field=id_tsk:integer' # IDX_ID
uri += '&field=name:string' # IDX_NAME
uri += '&index=yes'

layer = core.QgsVectorLayer(uri, 'name', "memory")

core.QgsMapLayerRegistry.instance().addMapLayer(layer)

When I look to the capabilities of the dataProvider (in QGIS Python console) I've got this:

>>> layer.dataProvider().capabilitiesString()
PyQt4.QtCore.QString(u'Add Features, Delete Features, Change Attribute Values, Add Attributes, Delete Attributes, Create Spatial Index, Fast Access to Features at ID, Change Geometries')

My data are contains in a QAbstractItemModel of my own.
So when the model data changed I want to update only element / geometry that has changed and not the whole layer (witch might be big).

self._layer.setReadOnly(False)
feat_2_del = set()
feat_2_add = []
feat_2_chg = {}
# +1 because we want the bottomRight row to be selected
for row_num in xrange(topLeft.row(), bottomRight.row() + 1):
    idx_geo = self._model.index(row_num, self._col_geo)
    wkt_v = idx_geo.data(QtCore.Qt.DisplayRole)

    obj = self._get_object_from_idx(idx_geo)
    feat = self._object_map.get(obj.id)

    geom = None
    if not wkt_v.isNull():
        wkt = wkt_v.toString()
        if wkt.size() > 0:
            # some object do not have any geometry defined
            geom = core.QgsGeometry.fromWkt(wkt)

    old_exist = feat is not None
    new_has_geo = geom is not None

    if old_exist and new_has_geo:
        # we just change the old feature
        feat_2_chg[feat.id()] = geom
    elif old_exist and not new_has_geo:
        # we remove the old feature
        feat_2_del.add(feat.id())
        del self._feature_map[feat.id()]
    elif not old_exist and new_has_geo:
        # we create a new feature
        feat = core.QgsFeature(id(obj))
        feat.addAttribute(IDX_ID, QtCore.QVariant(obj.id))
        feat.addAttribute(IDX_NAME, QtCore.QVariant(obj.name))
        feat.setGeometry(geom)

        self._object_map[obj.id] = feat
        self._feature_map[feat.id()] = obj
        feat_2_add.append(feat)

self._layer.dataProvider().deleteFeatures(feat_2_del)
self._layer.dataProvider().addFeatures(feat_2_add)
self._layer.dataProvider().changeGeometryValues(feat_2_chg)
self._layer.commitChanges()
self._layer.setReadOnly(True)
self._iface.mapCanvas().refresh()

When I add a geometry information to an object (addFeatures), the layer is correctly updated (with the according geometry). But If I change the geometry (changeGeometryValues) or remove a feature (deleteFeatures) nothing is done. But my dataProvider should handle it.

What am I doing wrong?


texts:
I use qgis 2.18. I've merged two lines into one with the merging tool from advanced toolset.
When I cut off a small part of the line, I see the original lines again?
In these two lines I have the same ID, and it should not be.
I try to make road network with start and end of roads nodes. In that reason Aeuquilibrae plugin cannot create Nodes on that merged line. This is not only one line, this is used like example.This is only the case with such lines. My question is: how to merge a line from several lines without ghosts?

The yellow line was the original line.
Picture below


texts: Raster calculator: using coefficients stored in csv table

texts:
I need to calculate an expression over 25 different rasters. The function is the same for all rasters, but the coefficients are raster-specific. I would like to tell PyQGIS to select the relevant coefficients and then feed them into GDAL raster calculator formula.
To be more specific, I have rasters named as nl_y, where y is a numeric identifier, and coefficients relative to each y stored in a table (my_file_csv).

What I have tried to do - following Houska's advise - is to: (i) import the csv table as a layer without geometry; (ii) create a string reporting the expression to be applied in the raster calculator; (iii) feed this string into GDAL raster calculator.
Below is the code doing this:

## (i) Import csv with coefficients
uri = "file:///my_file_csv?delimiter=,"
L_coeff = QgsVectorLayer(uri, "csvLayer", "delimitedtext")
QgsProject.instance().addMapLayer(L_coeff)

for y in year:
    ## (ii) Select coefficients and compute formula
    for f in L_coeff.getFeatures():
        if f["Var name"] == "Var 0":
            c0_y = f["Column " + y]
        if f["Var name"] == "Var 1":
            c1_y = f["Column " + y]
        if f["Var name"] == "Var 2":
            c2_y = f["Column " + y]

    formula_y = "(" + str(c0_y) + ")* (A > 0) + (" + str(c1_y) + ")* A + (" + str(c2_y) + ")* (A)^2"

    ## (iii) Raster calculation
    Fname_y = "nl_" + y + ".tif"
    processing.run("gdal:rastercalculator", {
        "INPUT_A": input_dir + "\\" + Fname_y,
        "BAND_A": 1,
        "FORMULA": formula_y,
        "OUTPUT": output_dir + "\\" + Fname_y})

I get output rasters but they have no-sense values - i.e. between 1.79769e+308 and -1.79769e+308, instead of having 0-63 range. My guess is that GDAL raster calculator does not correctly read the formula, even though when I print it, it is exactly identical to what I would write manually in GDAL algorithm.

Any idea on how I can apply raster calculation (not necessarily on GDAL!) in an automated manner using Python console in QGIS?


texts: ST_Latitude and ST_Longitude functions do not return exact values in MySQL

texts:
I tried to collaborate on this question from the Stack Overflow site.

The ST_Latitude and ST_Longitude functions (equal ST_X and ST_Y) of MySQL (version 8.0.17) returns a non-exact value for certain coordinates (e.g.: -31, -61 or -62). At first I thought it could be a bug in this MySQL functions and I reported it in bugs.mysql.com, however, from MySQL they indicate that it is not a bug and that I must read in depth the GIS standard.

DROP TABLE IF EXISTS `tbl_table`;

CREATE TABLE IF NOT EXISTS `tbl_table` (
  `random_point` GEOMETRY
);

SET @`SRID` := 4326;

INSERT INTO `tbl_table`
VALUES
  (ST_GeomFromText('POINT(-30 -31)', @`SRID`)),
  (ST_GeomFromText('POINT(-32 -33)', @`SRID`)),
  (ST_GeomFromText('POINT(-60 -61)', @`SRID`)),
  (ST_GeomFromText('POINT(-62 -63)', @`SRID`));

SELECT
  ST_AsGeoJSON(`random_point`) `geojson`,
  ST_Latitude(`random_point`) `latitude`,
  ST_Longitude(`random_point`) `longitude`,
  ST_X(`random_point`) `x`,
  ST_Y(`random_point`) `y`
FROM `tbl_table`;

geojson                                          |           latitude |           longitude |                  x |                   y
:----------------------------------------------- | -----------------: | ------------------: | -----------------: | ------------------:
{"type": "Point", "coordinates": [-31.0, -30.0]} |                -30 | -30.999999999999996 |                -30 | -30.999999999999996
{"type": "Point", "coordinates": [-33.0, -32.0]} |                -32 |                 -33 |                -32 |                 -33
{"type": "Point", "coordinates": [-61.0, -60.0]} |                -60 |  -61.00000000000001 |                -60 |  -61.00000000000001
{"type": "Point", "coordinates": [-63.0, -62.0]} | -61.99999999999999 |                 -63 | -61.99999999999999 |                 -63

See dbfiddle.

I understand that MySQL, MariaDB and PostgreSQL use the Open Geospatial Consortium (OGC) specification and after reading the specification anyway I cannot determine the reason why some coordinates do not return the exact value in MySQL.

Both MariaDB and PostgreSQL return the exact value of the same coordinates as demonstrated in the following examples: MariaDB (version 10.4.8) and PostgreSQL (version 11.5).

Has anyone faced the same problem before?.


texts:
I have a TIFF image with no coordinate reference system and I need that image to be in EPSG 32631. Is there a way to do this using GDAL?

Note: I am very new to Python


texts: Adding control to wfs layer disable panning and zooming on double click

texts:
I can zoom and pan wfs layer when I don't have added control but I can't pan the map or zoom the layer with double click when I add a control to map: code for control is :

var select = new OpenLayers.Control.SelectFeature(district);
    map.addControl(select);
    select.activate();
    district.events.on({
        featureselected: function(event) {
            feature = event.feature;
            feature.popup = new OpenLayers.Popup.Anchored
                    ("dis_pop",
                            feature.geometry.getBounds().getCenterLonLat(),
                            new OpenLayers.Size(285, 250),
                            '

' +
                            '
' +
                            '
' +
                            '' +
                            'Item' +
                            '	' +
                            '' +
                            'Value' +
                            '	'
                            '
' +
                            '

',
                            null,
                            true
                            );
            map.addPopup(feature.popup);
        },
        featureunselected: function(event) {
            var feature = event.feature;
            map.removePopup(feature.popup);
            feature.popup.destroy();
            feature.popup = null;
        }
    });

The code is working file but I can't zoom on double click. I also can't pan the map with mouse cursor.


texts: I would like to add a graphic/image on the top left corner of my view, would it be possible to make it static, meaning that wherever we zoom to or pan, the graphic will still stay put.


texts:
I am unable to edit WFS layers in QGIS. I have set up a service access rule in Geoserver to allow users associated with role "staff_edit" to edit WFS layers. I followed this tutorial to ensure that I did not miss a critical step. I am currently using QGIS 2.16 and Geoserver 2.9.0.

I have the following service rule in place:

wfs.*   GROUP_ADMIN,ADMIN,staff_edit

No other rules have been specified (yet). I am able to connect to the server and add WFS and WMS layers in QGIS, but I am unable to edit WFS layers. The following have been verified:

WFS service level set to "Complete"
"staff_edit" role given WRITE access mode
The data come from a PostGIS data store and credentials grant editing privileges (currently using admin credentials).
tried uploading a vector dataset that was not in PostGIS, and same inability to edit.

How can I allow a specified role to edit WFS layers in QGIS?

UPDATE: The problem only exists on QGIS 2.16.3. Editing works on QGIS 2.14.


texts: using arcpy to convert labels to annotation

texts:
I am trying to use ArcPy to convert a specific shapefile in my TOC to annotation in the map. I am not so much in the need for the tiled annotation tool that I've seen talked about as this is going to be used at the layout view at the same scale always. Essentially I am trying to label all the areas in the [area] field and then convert to annotation right after so they can be dragged around for best look as there can be sometimes up to 50 different labels.

I've been trying this with no luck

def onGenerate(self, event):
    dlg = wx.MessageDialog(self, 'labels generated','success',wx.OK | wx.ICON_INFORMATION)
    dlg.ShowModal()
    dlg.Destroy()
    mxd = arcpy.mapping.MapDocument("CURRENT")
    layer = arcpy.mapping.ListLayers(mxd, "project_app")[0]
    if layer.supports("LABELCLASSES"):
        for lblclass in layer.labelClasses:
            lblclass.showClassLabels = True
            arcpy.ConvertLabelsToAnnotation_defense(mxd, layer)
    lblclass.expression = " [NAME]"
    layer.showLabels = True
    arcpy.RefreshActiveView()


texts: Given a lat/long how to find the nearest node in a route calculated by pgr_dijkstra?

texts:
I'm creating a table route_<route_id> for the route given by pgr_dijkstra, like this. In my application, I have a websocket with a mobile client that it's sending it's current position every x seconds. How can I find the nearest node from the table route_<route_id> and return a subset of the route, from the current position to the destination?

The SQL I'm using to generate the route table is:

CREATE TABLE route_route_id as SELECT seq, cost,
             length::double precision, geom,
             ST_AsGeoJSON(ST_LineMerge(geom))::json->'coordinates' as coordinates,
             ST_AsGeoJSON(geom)
             FROM pgr_dijkstra('
                SELECT gid AS id, source::integer,
                       target::integer, length AS cost
                           FROM "ways",
                               convert_to_lon_lat(lat, long, '"ways_vertices_pgr"'),
                               convert_to_lon_lat(lat, long, '"ways_vertices_pgr"'),
                               false, false) a LEFT JOIN ways b ON (a.id2 = b.gid);


texts: How to combine two different raster attribute tables into one table? There is no common value.


texts:
I have three spatial datasets - mitch_ibra has about 2700 polygons, hev has ~1,200,000 polygons and prop_p1 with about 180,000. For each polygon in mitch_ibra I want to know how many polygons of prop_p1 are inside, with a minimum size of 5ha or overlaps the mitch_ibra polygon by 5ha. Finally I want the number of these polygons intersect polygons from the hev table.

WITH region AS
  (SELECT geom
   FROM mitch_ibra
   WHERE sub_name_7 = 'Tingha Plateau'
     AND lscape_nam = 'Dumaresq Channels')

SELECT COUNT(props.id)
FROM
  (SELECT DISTINCT initprop.id
   FROM
     (SELECT id,geom
      FROM prop_p1
      WHERE st_area(ST_transform(geom,3308))>50000) AS initprop
   JOIN region ON st_intersects(initprop.geom,region.geom),
     (SELECT hev.geom FROM hev JOIN region ON ST_Intersects(hev.geom,region.geom)) AS hevsel
   WHERE CASE
             WHEN ST_Within(initprop.geom,region.geom) THEN ST_Area(ST_Transform(initprop.geom,3308))
             ELSE ST_Area(ST_Transform(ST_Intersection(initprop.geom,region.geom),3308))
         END > 50000
         AND ST_Intersects(initprop.geom,hevsel.geom)) AS props

Above is my query so far, which gives me a processing time of about 658s for a region with 1300 polygons selected in the prop_p1 subquery. All have tables have GIST indexes.  How can I optimize my indexes/query to reduce processing time.  Is this as good as I'm going to get? (I have to run this over 300 different selections in the WITH statement)


texts:
i tried create simple code for searching in JSON point feature layer (us_states) but i dont know, why it is not working. What can be wrong with this code?

        // Map data and base layer
        var data = us_states;
        var map = new L.Map('map', {zoom: 5, center: new L.latLng([37.8, -96]) });

        map.addLayer(new
        L.TileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png'));  //base layer
        var markersLayer = new L.LayerGroup();
        var featuresLayer = new L.GeoJSON(data, {
                onEachFeature: function(feature, marker) {
                    marker.bindPopup('
'+ feature.properties.OBJECTID +'
');
                }
            });

        map.addLayer(featuresLayer);
        map.addLayer(markersLayer);
        var controlSearch = new L.Control.Search({
            position:'topright',
            layer: markersLayer,
            initial: false,
            zoom: 12,
            marker: false
        });

        map.addControl( controlSearch );
        for(d in featuresLayer) {
            var title = featuresLayer[d].OBJECTID,  //value searched
                loc = featuresLayer[d].loc,     //position found
                marker = new L.Marker(new L.latLng(loc), {OBJECTID: OOBJECTID} );//se property searched
            marker.bindPopup('OBJECTID: '+ OBJECTID );
            markersLayer.addLayer(marker);
        }


texts: Identify polygon by property attribute in Leaflet

texts:
On a Leaflet map, I need to get the center using polylabel of a polygon which needs to be identified by geo_id.

Codepen here.

A snippet of what I have so far:

$(document).ready(function() {
    $(".button").click(function(){

    var parcelNumber = $(this).attr('name');
    var matchingPolys = sections.filter(function (feature){
              return feature.properties.geo_id === parcelNumber;
            });
            var singlePoly = matchingPolys[0];
      var pcenter = polylabel(singlePoly, 1.0);
            console.log(parcelNumber);
            console.log(pcenter);
    });
});

This gives me an error of polygon is undefined from polylabel, which means the loop isn't selecting or outputting the polygon correctly.


texts: Python - Correct way of closing ogr.DataSource() objects

texts:
I learned that I need to close objects that I open, usually using a with open(): statement that does the closing for you after leaving it. With GDAL/OGR, this does not work and objects must be closed manually. In GDAL, this works by
from osgeo import gdal
ds = gdal.Open(MY_FILE, gdal.GA_ReadOnly)
ds = None

With OGR, I found two ways to do that:
from osgeo import ogr
ds = ogr.Open(MY_FILE)
ds = None

# outputs:
# R6025
# - pure virtual function call

# crashes the debugger

ds.Destroy()

# same

While it seems to work fine with both methods when I run the full script, it crashes my debugger when I run it in debug mode (PyCharm, in my case). In the console, ds = None seems to work, while ds.Destroy() crashes the whole console with the Pure virtual function call output.
Is this specific to PyCharm or am I doing something essentially wrong?


texts:
I am trying to write my first ArcPy script to iterate a Band Math process with requires a few steps. The uniterated, working script for a single group of bands in one Landsat scene is:

#import modules arcpy and string
import arcpy, string

#import the environment settings and spatial analyst extension
from arcpy import env
from arcpy.sa import *

#Check out the Spatial Analyst extension (must be available to check out!)
arcpy.CheckOutExtension('spatial')

#Add your workspace here
env.workspace = r'C:\WORKSPACE'
env.overwriteOutput = True

#Define inputs from TIFFs
Band2 = 'LANDSAT_B2.tif' #Band 2
Band4 = 'LANDSAT_B4.tif' #Band 4
Band5 = 'LANDSAT_B5.tif' #Band 5

#Band Ratios
R42 = arcpy.sa.Divide(Band4, Band2)
R52 = arcpy.sa.Divide(Band5, Band2)
R54 = arcpy.sa.Divide(Band5, Band4)
print 'Ratioed Bands'

#Slopes
a = 52
b = 11
c = 33.5
d = 120

#Algorithm processing
ALG = a - (b * R42) + (c * R52) - (d * R54)

#Saving final product
result = 'ALG.tif'
ALG.save(result)
print 'Success'

My goal is to incorporate these steps into a for loop, so that I can run a script on a folder or set of folders rather than retyping each different scene's filename into the beginning of the script.

My plan is to separate the original bands, which are individual rasters, into folders and call them as raster objects:

arcpy.env.workspace=r'C:/WORKSPACE/Band2'
B2rasters=arcpy.ListRasters('*.tif*')
arcpy.env.workspace=r'C:/WORKSPACE/Band4'
B4rasters=arcpy.ListRasters('*.tif*')

and then begin by iterating the band ratio process through the folders, so that the first band 2 in the band 2 folder is ratioed with the first band 4 in the band 4 folder, and the second band 2 and 4 are ratioed, etc:

for idx in range(0, len(B2rasters) - 1):
    b2 = B2rasters[idx]
    b4 = B4rasters[idx]
output = b4 / b2
result = output + idx
output.save(result)
print result

I don't feel confident in my ability to construct the iteration even after trawling through many examples of shorter algorithm iterations, as none of the examples I found use as many different rasters for multiple steps.


texts: Problem with GDAL distribution and netCDF support

texts:
I have a problem with a GDAL distribution downloaded from http://www.gisinternals.com. It is the "release-1800-x64-gdal-2-2-1-mapserver-7-0-6.zip", downloaded just today.

In the information about this distribution it is indicated that it includes support for netCDF format.
But when I run "gdalinfo --formats", netCDF  is not in the list of the supported formats. Moreover, when trying to open a netCDF file from a C# program that uses the gdal/csharp libraries of this distribution, it says  that it doesn't recognizes the format. Same thing when using gdalinfo passing a netcdf file (checked as correct with another program) as parameter.

If this is a problem of the distribution, where can I find a distribution of GDAL binaries with full support (netCDF, HDF...) and with the libraries for C#?


texts: How can I create a Vector layer from selected features with PyQGIS?

texts:
In my QGIS plugin I select dynamically features from a Vector layer.
and currently i create a new layer (Shapefile) from selection on combining all feature into a new feature:

theField = QgsField
self.theString = (os.path.expanduser("~")+'\.qgis\\statsRectangle.shp')
feat = QgsFeature()
shapeLayer2 = QgsVectorLayer(self.theString, "Upstream Area Of Interest", 'ogr')
geomtotSubwatershed = QgsGeometry.fromWkt('GEOMETRYCOLLECTION EMPTY')
nodLayer.setSelectedFeatures(selectFeatureIDlist)
UpstreamGeometry = QgsGeometry.fromWkt('GEOMETRYCOLLECTION EMPTY')
for elem in nodLayer.selectedFeatures():
   UpstreamGeometry= UpstreamGeometry.combine(elem.geometry())
feat.setGeometry(UpstreamGeometry)
QgsMapLayerRegistry.instance().addMapLayer(shapeLayer2)

but combining selected feature take a long time and froze QGIS UI.


texts: If I create Iso-Area polygons using QNEAT3 plugin in QGIS, is it possible to get an output and read it in R for doing further data analysis? If yes, how can I do it?


texts:
I am trying to add data into open jump from a table in pgadmin. However there no geometry column in the table . I do have geography column. Can I transfer the geography to geometry ( i believe it ST_Transfer() method ) . Then add to open jump.

Additional information. There are other table within the database but they already have geom column.


texts: I need to get a count of neighbors within 100m for each point in a set of 250 million points. The result would have the id of each point accompanied by the number of neighbors around them. Can this be done in postgis in a few hours (the table has an indexed geometry column)? How long should it take and what would be the most optimized query?


texts: Disabling ArcGIS Server web adf calls when zooming with a mouse wheel?

texts:
When zooming using the mouse wheel on a web adf based map, every time the map wheel clicks, requests are sent to ArcGIS Server for maps etc. (as evident in Fiddler).

Is there any way to harvest an event that would disable the requests until you got to zoom level you want?


texts:
If I use ArcGIS Desktop applications like ArcMap and ArcCatalog, at the ArcEditor/Standard, ArcView/Basic or ArcInfo/Advanced license level, is it using ArcGIS Engine under the hood?

I think the answer is no, but please enlighten me.


texts: Find out if a point exists in an irregular shaped ILayer

texts: Assuming I have an layer from a shapefile which is not shaped like a rectangle how would I determine if a point lies within this layer?


texts: How to set LocateMe Widget static or predefined location(or x & y)

texts:
I want to have locateMe button with static predefined x and y coordinates.

I have found this:

myWidget = new locateButton({
theme: "locateButton", // (optional). class name for the widget container. default: "locateButton"
map: myMap, // (required) map object. default: null.
visible: true, // (optional) show the widget. default: true.
highlightLocation: true, // (optional) show a graphic on the location when geolocated. default: true.
scale: null, // (optional) scale in meters to geolocate to. Defaults to accuracy from geolocation.
symbol: new PictureMarkerSymbol(require.toUrl("esri/dijit") + '/images/blue-dot.png', 21, 21), // (optional)
infoTemplate: null, // (optional) popup info template. default: null
useTracking: true, // (optional) use watchPosition instead of getCurrentLocation. default: false
setScale: true, // (optional) set scale on locate. default: true.
centerAt: true, // (optional) center at point on locate. default: true
geolocationOptions: // (optional). default: { maximumAge: 0, timeout: 15000, enableHighAccuracy: true}}, "locateButton");myWidget.startup();

But I can't make it run.

Any help how to achieve the LocateMe button with predefined values?


texts: Manually triggering delete using Sketch Widget for ArcGIS JavaScript API

texts:
I am using the sketch widget for ArcGIS JavaScript: https://developers.arcgis.com/javascript/latest/api-reference/esri-widgets-Sketch.html

There is a "delete" button after you have selected one or more graphics. Is there a way to trigger the delete operation through the API?

I don't see one in the documentation.


texts: Workflow copy attribute in Geocortex Viewer

texts:
Is there a way to copy a attribute from a feature and paste it into a other feature attribute with a workflow?

ArcMap has a internal Tool spatial adjustment where this idea is possible with the button Attribute Transfer tool.

Is this possible in geocortex viewer with a workflow?


texts: Layer Tree Collapse Alongside OpenLayers3

texts:
I have the following JSFiddle that includes a layer tree for my OpenLayers3 map.

It is currently set up so that when you click on the branch 'parent', the layers collapse. I'm looking to make it so that on the page load, the "Cities" and Background Mapping" 'parents' are collapsed and the user can choose which branch they expand.

Preferably I wouldn't like it to collapse all the way up to "All Layers". As per the comments below, Twitter Bootstrap may be used.


texts: qgis2web categorized renderer symbology with svg markers issue

texts:
In a simple QGIS (2.14.5LTR) project I've set up categorized renderer symbology for a small shapefile which looks like this:

Categorization is based on a text field, cp. my comment on qgis2web export categorized renderer symbology

Exporting this to leaflet with qgis2web (1.19.0) results in

In my opinion this is not an svg issue, because with single symbol renderers the layer symbology looks absolutely fine in the leaflet project, as proven by the screenshot below:

svg symbols are from https://openclipart.org/

The svg files in this case reside in the same folder as the shapefile and the qgis project.

The svg files, shapefile and the project can be obtained here: https://1drv.ms/f/s!AoRm0zIUU8Gdp1e975aOcc9Ksrkx

What could be my problem?


texts: get lat long from maps and displayed it to text box in java script

texts:
i make google maps autocomplete address textbox to get latitude and longitude without map by using google.maps.places.Autocomplete class we can implement Google places autocomplete textbox and display selected address latitude and longitude details

I've tried like this

                                location:<input type="text" id="txtautocomplete" style="width:200px" placeholder="enter the adress"/>

but it didn't work for me, the autocomplete address didn't show please help me


texts: Failed to load the GDAL native libs for GeoServer 2.15.0

texts:
2 days searching and looking for a guideline to enable GDAL extension for GeoServer, I found many questions here and I followed all steps, but unfortunately the problem still exist the default raster extension not changed.

I just need a clear guideline how to enable GDAL, ECW, MrSID extensions for GeoServer 2.15.0

edit:
error message

WARNING: Failed to load the GDAL native libs. This is not a problem unless you need to use the GDAL plugins: they won't be enabled.


texts: GetFeature EqualTo filter incluiding substring function?

texts:
I have the following table as a layer in geoserver:

I'm trying to get the features with the first 3 characters of the codigo field matching the Literal.

When doing a GetFeature request with EqualTo filter I get features even if the 3 first codigo characters are not equal to the literal. Note that codigo is an string. This is my filter:

<Filter xmlns="http://www.opengis.net/ogc">
    <PropertyIsEqualTo>
        <Function name="strSubstring">
            <PropertyName>codigo</PropertyName>
            <Literal>0</Literal>
            <Literal>3</Literal>
        </Function>
        <Literal>201</Literal>
    </PropertyIsEqualTo>
</Filter>

Is there something wrong with the filter?


texts: Skyline Tool: Using angle of elevation from observer point?

texts:
Situation:

I want to determine if I can aim an antenna at a satellite from a given location and still receive satellite signals without surface obstacles (e.g., mountains) blocking the antenna. The antenna azimuth would be 146° (magnetic) and the antenna elevation would be 19°.

Question:

I believe the Skyline Tool holds the answer, except that I do not know how to account for the 19° elevation of the antenna within the tool. I've attached a sketch of what I think may hold the answer: taking the sin of the antenna angle and the distance to the nearest mountain, determine the elevation in meters that the antenna must clear and compare this to the height of the peaks. The problem is that foreground foothills and lesser peaks may still obstruct the angle of the antenna, and I want to conduct this exercise on more than one peak, but all surrounding surfaces from 140° to 150° (magnetic).

This question 1by1 was asked at geonet.esri, but I have got exactly same problem.


texts: Changing string to integer using ArcGIS Field Calculator?

texts:
I have a table with two fields. One field (type short) is empty. One field (type string) is full of digits. There are no random characters in that field, just digits. I want to convert the string digits into integers and put them in the empty field.

How do I do this in the field calculator?

Figured it out. The strings of digits were 8 digits long. That would make an integer too large for a type short. I made a new field of type long and field calculator converted it easily.


texts:
I have a shapefile containing a lot of lines which I'm trying to separate based on direction.

Is there any way to calculate the line bearing against a compass (in degrees or just as North / North-East / East etc.) add this calculated value as a column to the shapefile?


texts: GeoNode commands: preceded by "geonode" or "python manage.py"?

texts:
It´s not really a problem but I would like to know the difference. I´m not a programmer. When I started with GeoNode 2.4 I used to run the commands as  geonode name-of-the-command
Migrating to 2.6, some of the commands didn´t work that way and now I´m using them as
python manage.py name-of-the-command

If I run geonode updatelayers, I get errors; but if I run python manage.py updatelayers, everything goes smooth.
Why?


texts: how to transfer a geotif with WGS84 to a mercator array for ploting it with matplotlib.basemap

texts:
I want to plot data from a sensor raster file (tif) to a map using basemap.
I have a geotif with 'WGS 84'

>>> gtif = gdal.Open("example.tif")
>>> print gtif.GetProjectionRef()
GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],

My map is using 'merc'

>>> m = Basemap(projection='merc', llcrnrlat=-26, urcrnrlat=48, llcrnrlon=-28, urcrnrlon=42, lat_ts=20, resolution='c')

How can I convert the tif in the merc projection?


texts:
I am making an app to store GPS and mobile sensor data.
My problem is how to store the data. The accelerometer data has a rapid frequency and the GPS data comes irregularly. Right now I've been storing it as a CSV file, but with this logic I get many redundant positions, because the sensor data is stored with a higher frequency

    id lat lon speed x y z time
    1  59  10  5     9 1 2 xxx
    2  59  10  5     5 4 3 xxx

The file will look something like this.

Is it possible to use a database logic with for example Postgis or do anyone have a good solution for this?


texts: Adding points to input point features option using spline method in ArcGIS Desktop?

texts:
I want to interpolate raster surface from points. I have excel file consisting of latitude, longitude and values( which I want to interpolate).

How do I add these point values in input point feature to extract the desired result?

Can I directly add points from excel sheet?


texts: How to regrid a raster after replacing the values of raster by the mean of 2 by 2 window and move to the next grid using R?

texts:
I am working on a raster image which is of 120m resolution. I need to make the raster divided into grids or windows of size 2 by 2 and then take the mean of the values in the grid and replace the cells by their mean by which the resolution will change to 240m (2 by 2). There should not be overlapping in the grids or window. The extremes can be padded to get the mean if the raster is having an odd number of rows or columns. The NaN values can be ignored and the results should be rasters.

How can I achieve this using R? I am quite new to it. I have converted it into a matrix and run a loop through it, Focal and sliding window and rollApply do not solve my problem.


texts: How to interpolate a point vector layer inside a boundary (polygon layer) QGIS/GRASS/SAGA?

texts:
I have to create a heatmap of water temperatures inside a lake for which I have sample points. The map should look much like the climate maps generated by climate change models. I have the temperature points and they were attached as attributes to a lat/long point shapefile layer.

I would like to interpolate these points using temperature as an interpolation attribute, but I'd like to contain the values within the lake.

I would prefer not to do an interpolation with the lake area polygon as a masking/clipping layer, but rather the interpolation should take in account the boundaries much in the way ANUDEM seems to do it for hydrological maps, if my understanding is correct. That is borders should influence the value of the interpolated value.

Please let me know if such interpolation is possible in one of: QGIS, GRASS, SAGA.

If ArcGIS can provide such functionality while the others don't I can get a temporary license to do this.

Edit: After doing more research I think that what I need is called interpolation with fault lines or barriers. I found an Open Source program surfit, which apparently does that. I am curios if in any of the software mentioned above there is something that does that?


texts: How do I get all the Layer Vectors added to a Map in OpenLayers 3?

texts:
I had a look at this post: How to get all Vector layers from OpenLayers map?

The following construct is shown:

    map.getLayersByClass("OpenLayers.Layer.Vector")

Unfortunately getLayersByClass is not recognised in this version of OpenLayers.

I add each of my vector layers here:

    if (!group.olSourceVector)
            group.olSourceVector = new ol.source.Vector({});

        if (!group.olLayerVector)
            group.olLayerVector = new ol.layer.Vector({ source: group.olSourceVector });

    this.olMap.addLayer(group.olLayerVector);

Then features are added to each LayerVector.

Just prior to drawing the map, I want to run through all the Layer Vectors to calculate the extent of the map.


texts: Identify Results tool not previewing fid instead of 'Display Field' in some cases

texts:
When I click on the "Identify Results," some layers of my GeoPackage display the attribute selected under layer properties/Display, but other layers will show the fid instead of the field I've attribute I've selected under layer properties/Display. Any idea why that might be?

Are there certain characters that could throw it off? Or maybe if there is another layer that has a feature with the same value in their respective display attributes?

In the below, the (53) shown from Riser layer is the fid, and not the "Name", but the others are all working as expected.


texts: Excel worksheet for Coordinate transformation between EPSG:31467 and EPSG:25832?

texts: for offline use I have been asked to provide a coordinate transformation worksheet in excel that converts between EPSG:31467 and EPSG:25832 and also WGS84. Has anyone any experience with something like this?