texts: Extracting Survey123 data from SQLite database

texts: I have a Survey123 survey SQLite database I need to extract data from. The original survey file accidentally got deleted and I'm trying to use the ReadS123DB tool to extract the data. However, whenever I try and run the tool either through the CMD or from Python IDLE the tool acts like it runs, but no output is generated. I've tried this with multiple datasets and manually specifying the input, but to no avail. Not sure if this is related to my own python incompetence or some update which invalidated some line of code in the tool. I've also tried Importing Survey123 SQLite Database from Device, but since my survey contains repeats and nested objects I don't think it would work regardless.


texts:
I aim to make the fusion of Landsat and Modis using STARFM at USDA StarFM
I tried to conduct this process in R, in which we call the software using this code:

system2(command="./StarFM.exe",args="StarFM_config.txt", wait=TRUE)

However, this cannot run to produce the final results.
I also tried to run the software from Command Prompt as shown below:

cd .\StarFM\source
StarFM.exe

And got this result:

Can you help me to check, and run this software in R, especially to process Remote sensing data?


texts: Easiest way to extract all city|town|village out of OSM?

texts:
I'm not a techie and I don't have the space to download the 500GB planet OSM file and run through osmosis.  What is the easiest way for me to download all nodes of type City+Town+Village (as defined by OSM), pulling just these parameters:
•    Node ID
•    Node Latitude
•    Node Longitude
•    Key (city or town or village)
•    Name
•    Population
•    is_in:state_code
•    is_in:country_code
Please help (in a way a non-techie can understand).  I've come across literature on many ways to pull data from OSM but instructions haven't been clean enough for me to figure out.


texts:
I'm using pyqgis within an standalone application.
Now I wanted to create a custom expression to to use as filterExpression in rule of ruleBasedRendererV2. (I tried the renderer with existing functions before and everything worked fine!)

Custom qgisfunction (bitwise AND-function - as there is not such a function yet) looks like this:

from qgis.utils import qgsfunction

@qgsfunction(args="auto", group='Python')
def bitwise_and(value1, value2, feature, parent):
    return value1 & value2

I also tried to register it, via:

from qgis.core import QgsExpression
QgsExpression.registerFunction(bitwise_and)

Testing the expression via:

expressionString = 'bitwise_and(2, 2) = 2'
exp = QgsExpression(expressionString)
if exp.hasParserError():
    print(exp.parserErrorString())

(Giving 2=2 so always True) works after registering function, and fails if not. But still, I can not use the expression within a rule...

Right now, the file containing the function lies on top in my project folder, which is automatically included to PATH, isn't it?! (According to Nathan Woodrow, file can be placed anywhere in PATH - see https://nathanw.net/2012/11/10/user-defined-expression-functions-for-qgis/).

UPDATE: Registering of function works so far, with code from above - now the problem lies in function itself...

UPDATE #2:

Very strange!
Simple function do_nothing getting two arguments works fine in my code:

@qgsfunction(args="auto", group='Python')
def do_nothing(value1, value2, feature, parent):
    return value1

Using expression expString = 'do_nothing(2, 2) = 2'

UPDATE #3:

Simple change of do_nothing to add_one does not work any longer:

@qgsfunction(args="auto", group='Python')
def add_one(value1, value2, feature, parent):
    return value1 + 1

calling it via expString = 'add_one(2, 2) = 3'.
That's why I think it has something to do with returning type?!


texts: Obtaining the portion of upstream cell from taudem

texts:
In taudem I know that I can use StreamNet function:

StreamNet -fel <felfile> -p <pfile> -ad8 <ad8file> -src <srcfile> -ord <ordfile> -tree <treefile> -coord<coordfile> -net <netfile> [ -netlyr netlayername] [ -o <outletfile>] [ -lyrname <layer name>] [ -lyrno <layer number>] -w <wfile> [ -sw]

pfile: Input D8 flow directions grid 'p'

srcfile: Input stream raster grid 'src'

ordfile: Output grid of channel network Strahler order 'ord'

ad8file: Input grid with D8 contributing area 'ad8'

felfile: Input pit filled elevation data grid 'fel'

treefile: Output textfile with list of links in channel network tree

coordfile: Output textfile with list of coordinates in channel network tree

outletfile: input outlets file (OGR readable dataset).This must contain a field named "id" that is

used to identify the links upstream of outlet points

layer name: OGR layer name if outlets are not the first layer in outletfile (optional)

layer number: OGR layer number if outlets are not the first layer in outletfile (optional)

Layer name and layer number should not both be specified.

wfile: Output grid of watershed identifiers 'w'

netfile: Output file (OGR writeable data set) of resultant channel network.

netlayername: layer name of netfile (Optional)

to obtain the streamline shape, by reading off the values in the output file treefile and coordfile.

But what about obtaining, for each cell, the portion of upstream cells that contribute flow to it?

I suspect that the answer lies somewhere in the command argument wfile( with the file extension .tif) and netfile ( with the extension .shp). But these are binary files, and I don't know how can I extract the information out.

Any idea?

Note: I am aware that shape and tif files can be opened by GIS suite, but that is not my purpose. I want to extract those information out because I want to be able to render them on my own choice of viewer, and maybe do some statistical analysis on it using statistical tools of my choice


texts: ArcMap, Joining information from lyr with attribute table

texts:
I'm working on a landcover project with ArcMap 10.1 right now and use the ModelBuilder a lot.

So far I extracted a table from a certain layer using the Summary Statistics tool.
However, one field of the mentioned table does only contain a code that isn't of much use to me. I do have a lyr file containing the right legend, but I cannot simply join the lyr file with my extracted table.

Is there some way to connect these sources of information? I do need to write some documentation on what tools I used during the project (using the ModelBuilder) and therefore don't want to resort to just write the correct code down manually within my extracted table.

Any help is highly appreciated!

Best Regards, T.


texts: Removing  From all Fields and all Feature Classes in a GDB

texts:
Have a gdb with 13 FCs and 70 Fields in each FC (same fields). Some fields have 'Null with <> around', which I would like to remove. (Not the field itself or data point. As in set null to empty "")

The only way I can think of doing this is by:

and supplycode provided by here by esri

However I would need to do this 70! times. Is there a more efficient way of doing it?


texts: Why does aggregate polygons consistently error out?

texts:
I'm using aggregate polygons on a dataset with three large study areas. It did not work. So I split it into separate areas and did it again. Two of the three areas (which were smaller than the third) aggregated correctly. The third I split into municipalities (since it had a few in it) and then aggregated each separately, apart from two, which refused to aggregate. The error message was not forthcoming, it just says "Background processing has encountered a serious error and is unable to continue. And then it allows me to send an error report to ESRI, which I have, though it's been almost innumerous times.

So I used grid index features and made a grid over one municipality and then split the data by that grid. I ended up having to do 5 x 5 km grids (making 42 index polygons) because anything bigger would not process properly. That actually sort of worked, except for five of the grids which would not. I tried vainly to do 'multipart to singlepart' which worked for some of them, though in this case, I had to further subdivide each 5x5 grid into 4 more indices. In that case, quite a few actually processed, but there is one or two stubborn grids which will not process at all and produce the same error.

I even went into each grid's data and used the editing toolbar to merge some of the polygons together manually, before using aggregate polygons.

Has someone else had similar trouble? Are there any workarounds or suggestions? Because I will otherwise have to repeat the process two or three more times (as it is a time-series project).

edit: I tried dissolve and it does sort-of do what we want, though is there any way to add a distance value to it like aggregate polygons has?


texts: How do I find out when a shapefile was created?

texts:
I am working on a project at an agency with a large amount of prior data stored in a  shared database connection .sde folder to which I have access. While this is incredibly helpful, I would like to know when these files were created so I can tell how recent/up to date the data is, in situations when the date is not noted in the properties > description for a file.

Is there a way to either (a) find out when a shapefile or database was created or (b) view a database connection folder in windows explorer?


texts:
I am passing in a geofeatureLayer to a method that creates a UniqueValueRenderer.
I am using the following code to sort that renderer so my values will be in some kind of order. This seems very slow, however. Is there a faster way to sort a featureClass?

IFeatureClass featureClass = geoFeatureLayer.FeatureClass;
ITableSort tableSort = new TableSortClass();
tableSort.Table = featureClass as ITable;
tableSort.Cursor = featureCursor as ICursor;
tableSort.Fields = classifyParms.ClassifyField.Name;
tableSort.set_Ascending(classifyParms.ClassifyField.Name, true);
tableSort.Sort(null);
featureCursor = tableSort.Rows as IFeatureCursor;


texts: Have Feature Click Events for Two Layers On, but only one is active at a time in CARTO.js 4

texts:
I am currently developing an application using CARTO.js 4, see here. We have a slider where layers for different years of data are loaded on both sides so you can compare them. We have a feature event where when you click on an area, a pie chart pops up. To get this to work I had to use two clients, one for each side and essentially copy everything from one side to another; this included the function we wrote for our feature event which brings up the pie chart.

  dominanceLayer2010.on('featureClicked', (f) => clickedOnFeature(f));
  dominanceLayer2017.on('featureClicked', (f) => clickedOnFeature(f));

  //functions removes popups that are not pinned when a new child is open
  function clickedOnFeature(featureEvent) {
    if ($('#popUpHolder').children().length > 0) {
      var popUpChildren = $('#popUpHolder').children()

      // console.log(popUpChildren.length)
      var count = popUpChildren.length;
      while (i = count--) {

        var pin = $(popUpChildren[i - 1]).find("svg")
        // console.log(i)
        // console.log($(popUpChildren[i - 1]).find("svg"))
        if (!$(pin).hasClass('pinned')) {
          $(popUpChildren[i - 1]).remove();
        }
      }
      // $('#popUpHolder').empty(); //commenting this out creates multiple pie charts

    }
    var originInfoNames= Object.keys(originInfo);
    for (index in originInfoNames){
      popFactory.pieChartData[index].value= featureEvent.data[originInfo[originInfoNames[index]].tableInfo.alias]
    }

    var popCount = popFactory.newPopUp(featureEvent);
    $("#popUp" + popCount).css({
      "left": event.clientX + 5,
      "top": event.clientY + 5,
      "visibility": "visible"
    });

    popFactory.headerText= featureEvent.data.latino_pop;
    var pie = new d3pie("pieChart" + popCount, popFactory.pieConfig());

  };

It seems the layer that gets turned on last is the one that the feature event works for. For example, we have the 2017 layer loaded on the right if you click on an area here and bring up the chart and then slide so the 2010 layer is showing and click on the same area, it will give you the data for 2017.

Anyone know why this is happening?

Here is the repository in case you need to see more:https://github.com/SBUtltmedia/cartoDominance

I just looked at this again and the problem is that the slider doesn't clip the data, if you click something that is not being exposed the function still fires. I need to figure out a way for the client to know which layer to fire the function for, maybe cursor position. I'm going to start investigating this.


texts: How to calculate the distance between point and a polyline along specific direction?

texts:
I want to find the closest road for point star. By using the Near tool in ArcMap, I can get the closet distance a since the closest road is road 1. However, following the red arrow direction, the closest road is road 2.

I am just wondering does anybody know any function in arcpy or other libraries that can calculate distance b.


texts:
How do you do raster division (ArcPy raster calculation)? I want to divide the “.tif” in one folder by the “.tif” in the other folder.
This is my code：
import arcpy
from arcpy import env

from arcpy.sa import *

arcpy.env.workspace = r"C:\Users\giser\Desktop\test1"

raster_list1 = arcpy.ListRasters("*")

print (raster_list1)

arcpy.env.workspace = r"C:\Users\giser\Desktop\test2"

raster_list2 = arcpy.ListRasters("*")

print (raster_list2)

try:

    if arcpy.CheckExtension("Spatial") == "Available":

        arcpy.CheckOutExtension("Spatial")

        print ("Checked out \"Spatial\" Extension")

    else:

        raise LicenseError

except LicenseError:

    print "Spatial Analyst license is unavailable"

except:

    print arcpy.GetMessages(2)

for raster1 in raster_list:

    ras1 = raster1

for raster2 in raster_list2:

    ras2 = raster2

    outraster = Raster(ras1)/ Raster(ras2)

    outraster.save(r"C:\Users\yutao\Desktop\test1" + "\\1" + ras1)*

RuntimeError                              Traceback (most recent call last)

<ipython-input-12-36e5c3a3bfbf> in <module>()

      for raster2 in raster_list2:

          ras2 = raster2

--->      outraster = Raster(ras1)/Raster(ras2)

          outraster.save(r"C:\Users\giser\Desktop\test1" + "\\1" + ras1)

RuntimeError: ERROR 000732: Input Raster: Dataset 2001maxndvi.tif does not exist or is not supported


texts: What is the most user friendly OSM route planning API?

texts:
I ran a workshop last week on GIS for Transport Applications. In Nick Bearman's tutorial (available from the previous link) we used Google's routing API (via ggmap) and CycleStreets.net routing algorithm via my new package stplanr. Here's an illustration of what the latter can do in R - love the fact this is all based on open data and open software:

I'm looking to include additional OSM routing services - which would people recommend?


texts: GeoJSON format (Geography) for posting in Postman tin PostgreSQL

texts:
I am trying to do an update posting on Postman for the first using PostGIS / Geography file type and I don't know where I am going wrong.
I have included the model below of the table in pgamdin 4 / PostgreSQL that I need the posting to go to as well as the json format. I know about GeoJSONs but I'm not sure if these are just JSONs with coordinates. I am trying to pull object description details along with the co-ordinates so i can display them on the marker.
Am I going about this the right way?
My JSON posting that isn't working is below:
{
    "diveSchoolName": "Forever Dive",
    "diveSchoolLocation": "Madagascar",
    "diveSchoolRegionID": 5,
    "diveCertificatesOfferedID": 1,
    "diveSpotsOfferedID": 1,
    "diveSchoolGeo": {
        "type": "Feature",
        "geography": {
            "type": "Geography",
            "coordinates": [-18.7669, 46.8691]
        }
    }
    }

model
module.exports = (sequelize, Sequelize) => {

    const approvedDivingSchool = sequelize.define("approvedDivingSchools", {

        diveSchoolID: {
            type: Sequelize.INTEGER,
            primaryKey: true,
            autoIncrement: true,
            allowNull: false
        },
        diveSchoolName: {
            type: Sequelize.TEXT
        },
        diveSchoolLocation: {
            type: Sequelize.TEXT
        },
        diveSchoolRegionID: {
            type: Sequelize.INTEGER
        },
        diveCoursesOffered: {
            type: Sequelize.INTEGER
        },
        diveSpotsOfferedID: {
            type: Sequelize.INTEGER
        },
        diveSchoolGeo: {
            type: Sequelize.GEOGRAPHY
        }},
        {
            timestamps: false
        })

    return approvedDivingSchool;
};

If I am posting a GeoJSON would it be inserted into one field in my PostgreSQL table or if I am pulling data from a table with coordinates is it automatically a GeoJSON? I'd prefer not to be duplicating work if information is already included in the GeoJSON.


texts:
I'm currently learn to build urban flood model using FLO2D plugin on QGIS. It's running well when using sample data.

But when I try to build model with my own data it doesn't work. I'm stuck on detailed graphic process for flood. I even tried to make a model by slightly changing the boundary from sample data to looking for the error, but it's not work again like before.

I learn from this tutorial: FLO 2D QGIS Start a Project.
Can anyone give suggestion to me? Because in my case, there is no error notification appear.


texts: Evaluating string in Arcade 'if' statement?

texts:
I'm tasked with building an Arcade expression that evaluates whether a string field includes given characters, and returns values accordingly. Attribute UTILITY_GROUP may be E, EU, G, or GU. To cut down on the number of comparisons required, I had hoped to do something like (pseudocode):

var utility = $feature.UTILITY_GROUP
if ('E' in utility) {
  return 'electric'
}
if ('G' in utility) {
  return 'gas'
}

Is there a way to express "if string includes [character]" in Arcade? I didn't see anything documented in the logical operators.


texts: Calculating distance between consecutive rows in QGIS


texts: PostGIS - Nearest Neighbour Join to sum values

texts:
I have a table of polygons T1 and a second table of points P1

Each point has a value

I want to sum the values of P1 for each T1, with P1 being the closest set of points to T1

At the moment the only way I can see to do this is compute a nearest neighbour for each P1 from the T1 table. Then having done this compute the sum. It seems inefficient as P1 contains millions of rows compared to T1 which has a few hundred thousand.

Is there a better approach? It is critical that points in P1 are not allocated to multiple polygons in T1. This is current approach

SELECT distinct on(B.id) B.id,
   sum("count") OVER (PARTITION BY B.id) as point_count
FROM P1
INNER JOIN LATERAL (SELECT id FROM T1 ORDER BY P1.wkb_geometry <->
T1.wkb_geometry LIMIT 1) B
ON TRUE


texts: How to write a (memory) layer from QGIS to PostGIS?

texts:
I am looking for an "elegant" solution to write, from a python script, a layer from qgis and into a new table of a PostGIS database. I thougth that this might be solved by QgsWriteVectorLayer, but I could not find any hint that this is supported. I only found the solution that this can be done e.g. with Qtsql, but then one has to:

analyse the layer (field names, field datatypes, srid, ... )
create the table and
fill the table with date.

Question: Does pygis offers a "simple" solution to write a complete layer directly into a new PostGIS table?


texts: Getting started with 'shell scripting loops' in GRASS

texts:
I'd like to work through the MODIS-Aqua example here.

But typically the hardest bit is getting off the ground.  What specifically does it mean by 'shell scripting'?  is it referring to the command console?  Python console?

To get started I tried  the 'Decompress' example bzip2 -d A20023532002360.L3m_8D_SST_4.bz2 in both command and python consoles - neither worked - I expect I need to get the file path in there, but not sure of the syntax. The file currently resides in the Project>Mapset folder where I thought the default path would lead - is that the best place for it?

The next big hurdle will be creating the loop to process hundreds of hdf files, but i can't get off the mark until I work out some of the very basics.  Sorry, this is my first foray into GRASS scripting.  Any help would be most appreciated!


texts: QGIS Export file Automatic Geometry Type issue

texts:
It appears that the new version of QGIS has an issue whereby if you save a new file (using 'Export > Save Features as) and then select 'Automatic' for the geometry type, it does not migrate the geometry type correctly.
For instance, I am saving a Mapinfo TAB file of geometry type Polygon (multipart). I perform the above export to a new file of same type and select the automatic geometry type.
The new file that is created SAYS 'Polygon (Multipart)' - however, attempting to save any records into this file throws an error message
"Geometry type is invalid for this format"
The work around is when exporting the file, to select the actual geometry type in the drop down on the export menu. Again, the output file is exactly the same and geometry format says Polygon Multi-part, but this time, saving records works.
Does anyone else experience this with file formats other than Mapinfo?
QGIS v3.16


texts: ArcPy: Aggregate points using multiple aggregation distances and output names

texts:
I would like to be able to run a whole load of aggregate polygons with a variety of aggregation distances, using Arcpy. I have successfully done it using model builder, but specifically want to do it using ArcPy as I am comparing the environments for research as a part of degree.
The code I have below is a start, but it errors out and was hoping someone would be able to help!

import arcpy
from arcpy import env

env.workspace = "E:/temp/Data.gdb"

in_point_features = "E:/points.shp"
out_feature_class = "E:/points"+aggregation_distance ".shp"
aggregation_distance = ["100 meters", "200 Meters"]

arcpy.AggregatePoints_cartography(in_point_features, out_feature_class, aggregation_distance)


texts:
I guess no one has not asked this question yet.

Let's say if I want to store my CSV to GDB feature table then the code below works just fine

ogr2ogr -f "filegdb" test.gdb my.csv

and I thought just creating an empty GDB is

ogr2ogr -f "FileGDB" test.gdb

and obviously it does not work.

Will there be a way to create an empty GDB using ogr2ogr?


texts: Long term NDVI mean from different resolution satellite imageries

texts:
I have 10 years satellite imageries of different resolution (Landsat TM, LISS3, AWiFS and Landsat 8). My aim is to generate long term mean NDVI from that data.

Since the resolutions are not same, i couldnt able to calculate the long term NDVI mean. How can i reduce the resolution of images with respect to one reference image.


texts: Does MySQL only support bounding boxes calculations instead of the real polygon geometry?

texts:
Within our company we are deciding which DBMS we will use to store spatial data. IT wants MySQL with the Spatial extension because there is already knowledge about handling MySQL DB. We are doing quite simple calculations with the data. The only thing I am worried about is the polygon handling of MySQL.

Is it still true that MySQL only supports bounding box calculations instead of the real polygon geometry? Or is this problem gone with newer versions of MySQL?


texts:
I have two vector layer files, 'hausnummer_cropped5' and 'gebaudepoints_cropped5 - Kopie2'. Both have a point geometry and I am trying to shift the locations of the gebaudepoints_cropped5 points onto the corresponding hausnummer_cropped5 points (there is a one-to-one correlation between them) if certain conditions are met (their original locations are slightly different). There are 68 gebaudepoints and 62 hausnummers.

The code I have written is this:

import arcpy

gebaudepoints_no = 68
hausnummer_no = 62
count = 0

with arcpy.da.UpdateCursor('hausnummer_cropped5', ['HAUSNR','STRASSE','SHAPE@XY']) as cursor1:

    for row1 in cursor1:

      with arcpy.da.UpdateCursor('gebaudepoints_cropped5 - Kopie2', ['ZA3',     'ZA2','GN5','SHAPE@XY']) as cursor2:

         for row2 in cursor2:

                if (row2[0] == row1[0]) and (row2[1] == row1[1]): # and (row2[2] > 0):

                  X = row1[2][0]
                  Y = row1[2][1]
                  count += 1

                  cursor2.updateRow([row2[0], row2[1], row2[2], [(X),(Y)]])
                  #cursor2.deleteRow()

                  print count, gebaudepoints_no - count

There are 2 things I don't understand:

After the code has been run, the value of 'count' is 92 and not 68. That means the loop has run 92 times even though it should have run only 68 times.
Most of the gebaudepoints actually get placed on top of the hausnummer points, but around 5 to 6 don't, even though they meet the required conditions. However, if I use the cursor2.deleteRow() statement (which is written as a comment above), all the gebaudepoints get removed so that indicates that the code is actually going through those points as well, but it's just not shifting them.

I didn't have these two problems when I was using a different set of gebaudepoints and hausnummers so I'm wondering why this is occurring in this case.

I am using ArcGis 10.2.2


texts: Simulate random error on data with pyqgis

texts:
I am developing a small tool to complement GPS data with mutual distances, and I need to simulate some input data. So I thought I would introduce a normally distributed error on the (known) positions, to simulate the GPS measurements. To make things easier, my simulated GPS measurements will have independent errors on the two coordinates, both averaging to 3m.

I want to do this within QGIS, and my code looks like this:

# just the standard initialization
#
from qgis.core import QgsApplication, QgsFeature, QgsGeometry
from qgis.core import QgsMapLayerRegistry, QgsProject

from PyQt4.QtGui import QApplication
from PyQt4.QtCore import QFileInfo

app = QApplication([])
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()

project = QgsProject.instance()
project.read(QFileInfo('/home/mario/Documents/GIS/heliconiaBB.qgs'))

# import normal distribution from 'random'
#
from random import gauss

# all layers are in WGS 84, or I want the code to be independent from the
# layer coordinates reference system

# the points in my virtual reality
reality = QgsMapLayerRegistry.instance().mapLayersByName("points")[0]
# the points as measured by my virtual GPS machine, to be simulated here
gps_readings = QgsMapLayerRegistry.instance().mapLayersByName("gps")[0]

# the two layers have the same fields set, including 'code'
fields = reality.fields()

gps_readings.startEditing()
pr = gps_readings.dataProvider()

for feat in reality.getFeatures():
    dx = gauss(0.0, 3.0)
    dy = gauss(0.0, 3.0)

    # the following line is just pseudo-code
    fake_gps_reading = feat.geometry().asPoint() + (dx, dy)

    new_feat = QgsFeature(fields)
    new_feat['code'] = feat['code']
    new_feat.setGeometry(QgsGeometry.fromPoint(fake_gps_reading))
    pr.addFeatures( [ new_feat ] )

gps_readings.commitChanges()

I do not know how to add the offset (in metres) to the point location, which is in degrees. I would be quite happy with the converse functionality as offered by QgsDistanceArea, which gives me the distance in metres between points on an ellipsoid given their latitude and longitude.


texts: Unable to install QGIS Plugins on Windows

texts:
I've installed QGIS 3.16 Hannover on my laptop about three days ago (Windows 10 Home 1909) and am unable to install any plugins, nor get access to any of the core plugins save four.

When going to settings, I get the ever-annoying "The settings on this tab are only applicable for Python Plugins. No Python support detected, thus no settings available" message.
Unlike most questions about this subject I've seen on here, updating didn't work; it still shows the same error.


texts: GeoServer Timestamp/Timezone error convert

texts:
Date format in PostgreSQL, -> 2019-10-30 08:46:44

Created WFS, but got a problem.

GeoServer is converting the date. I think it somehow something with timezones.

Output WFS -> 2019-10-30T07:46:44Z

So one hour difference. Need timezone Amsterdam, not the standard.

How to solve this?


texts: ArcGIS Mosaic to New Raster() tool seems to interpolate data gaps on one computer/version but not another?

texts:
i have a group of raster files in .tif format that I want to mosaic together into one .tif.

I've noticed that I get data gaps along the seams when I do Mosaic to New Raster() on my home computer (ArcMap 10.1). Also, these data gaps are not necessarily where there is no data because when I display the original .tif files, there is coverage where the final mosaic-ed raster shows "no data" values. There should be some gaps because the original rasters do not completely overlap, but there is more coverage than what is shown here as gaps. Here is a screen shot of the mosaic-ed raster with the data gaps:

Then, when I do the same operation on the same rasters in the same order, but on a different computer (ArcMap 10), there are no gaps at all. Here is a screen shot of the mosiac-ed raster in the same region as above:

Ideally, I would like a final single raster with no data gaps so I would like to use the no gap raster, but is it reliable? I've read that Mosaic to New Raster does not interpolate the no data areas, so what could be going on? How do I know which final layer is reliable?


texts: Generate a viewshed and merge it with another image using GRASS/GDAL command line?

texts:
I've found a tutorial showing exactly what I want to do (merge viewshed and hill shade or map) but unfortunately it is for QGIS users and I need it in pure command line instructions because I'm in a headless environment.

Must be only GDAL/GRASS command line and use SRTM HGT as DEM source and hill shade source.

I don't know QGIS enough to think it is very compact and directed to expert users or if it is just full of blah blah blah because it have a lot of images and almost none real instructions.

Can someone translate this tutorial to pure command line GDAL/GRASS so I can achieve the same result like the pictures below?


texts: How can the ogr2ogr help to merge a shapefile?

texts:
I have over 30 shp files which is named "street.shp" into the sub-folder. There is over 30 sub-folders. In the sub-folder, there is another shp files which may be named a building.shp...etc. I only want to merge the street.shp files. The structure of all street.shp is exactly same!

I tried this script -

for /R %f in (*.shp) do ogr2ogr -update -append merge.shp %f -f “esri shapefile” -nln merge

How can the script find the shp file named "street.shp" only? Thanks a lot.


texts: google geocoding service combined with openlayers - not working

texts:
I'm trying to combine google's geocoding service with openlayers. So when a user types an address and hits enter, map zooms there and renders a simple point. I edited my code according to the previous link. Does not work. I'm using a vector layer as basic (comes from Geoserver/PostGIS), a google's layer as second basic. I also added a simple vector layer to render the adrress's point. Do  I have to bind google's layer and geocoding service somehow?

Here is my code.
Can you spot what's wrong?

var geocoder;
function codeAddress()
{alert ('inside function');//works
var address = document.getElementById("address").value;

//tell to user, he did not type anything
if (address==""){document.getElementById("geocodMsgG").innerHTML = "You gave no input";}
else{

alert(address);//works

geocoder.geocode( { 'address': address}, function(results, status)
  {alert(status);//not work

   if (status == google.maps.GeocoderStatus.OK) {
     alert('inside 2 function');//not work
     alert(results[0].geometry.location);//not work
     //create a geometry point
     pointgc = new OpenLayers.Geometry.Point(results[0].geometry.location);

     //create a feature point
     feature_point_gc = new OpenLayers.Feature.Vector(pointgc);

       //lfl contains bounds that i set earlier in my code
       if (lfl.containsLonLat(results[0].geometry.location)){
       //render the point on the layer named geoco
       geoco.addFeatures([feature_point_gc]);
       geoco.refresh({force:true});
       //zoom there
       map.setCenter(results[0].geometry.location);
      }
      //if point is not inside the bounds (lfl) then dont add - print msg to user
      else{document.getElementById("geocodMsgG").innerHTML = "out of bounds";}
  }

  //refers to the second if(geocoding's status)
  else {document.getElementById("geocodMsgG").innerHTML = "error";}

});//geocode
}//major else (if address empty or not)
}//function


texts: Exporting table to XYZ ASCII file via ArcPy?

texts:
I'm looking for a way to export an ArcGIS table (created with the Sample tool) to a text file via ArcPy.

I can do this in ArcGIS via the context menu by right-clicking the table, but have not found a way to script this.


texts:
I'm doing a Python simulation of a vehicle, driving along a specified path. I generate dummy GPS readings, that fetch the current position of the vehicle, with some randomized error added. After a track is complete, I obtain a set of collected GPS points, that looks like this:

Assuming that we have no knowledge on how the physical path looks like (for performance reasons, we do not want to have anything to do with the map), and that we only know these points, I would like to detect the turns of this track. I want to do this, so that later, I can group these points by linear segments, and then construct a convex hull for each of them separately, like this:

One idea that I already had, is, given a set of data points, construct a piecewise linear approximation of these points, and then extract the coordinates of the resulting line. However, it is still somehow problematic, because I would have to look for the closest existing point from the dataset, which seems to be time-consuming. The best outcome for me would be to identify right away the specific point (or points) at which the turn occurs. Every track has a list of its recorded points, and I have to cut it into pieces somehow, but to do that, I need to know at which points it turns.

more exact convex hull problem definition
Maybe this approach is wrong in general, maybe there is a much easier way, to find a convex hull without large empty areas for a group of points right away?

Additional information:
The points are NOT an ordinary point cloud, they are stored as a list, which means that they are ORDERED.


texts:
We have an adobe illustrator file with nice base map from ArcMap that has a layer that we made edits to.  It is our town local streets with some additions and some subtractions.

we export the streets from ArcMap to .ai and it does not line up, apparently above file was rotated a bit and stretched, transformed or skewed etc. I can get it close, but so not perfect at all.

Is there any hope to line up our modified .ai layer to the original layer in the above file so it lines up perfectly?

We want to do this because the first .ai file has nice base maps etc, and just the one layer is changed.  WE also don't have access to the original ArcMap.


texts:
I am trying to add a geojson layer with interactive functionality to my leaflet map. I am mostly using this (dated)blog: http://palewi.re/posts/2012/03/26/leaflet-recipe-hover-events-features-and-polygons/

My result to date:

http://185.19.29.113:8080/cresh_map/mapgeojson.html

it changes the style of a selected feature on mouseclick
a popup come up on mouseclick
it changes the style back on mouseover

Could anybody help me with this question:

How can I change back the style on the next mouseclick instead of mouseover?


texts: Extra row and column in a raster after GDAL grid import

texts:
I'm exporting an ESRI grid to TIFF format and after that I'm reading this exported raster file into SAGA GIS. However when Saga reads this file in there appears to be an additional row and column in the same file and thus changed raster extent. Also cellsize has got some pseudo decimals (2 m -> 2.0000000000000022 m).

ArcMap (10.1) reports 24000 and 21000
while Saga (2.1.0) reports 24001 and 21001 columns and rows for the same file.

Coordinate system is EUREF FIN TM35.

Update:
Saga uses GDAL import grid -module as a default to reading rasters in. If one selects this option from module tree (not by load -> file) there are two options to select:

  select from multiple bands [ ]
  transformation [ ]

Under transformation there are different interpolation options, which seem to cause the problem. If none is selected, Saga reads my raster in with correct dimensions, but changes cellsize to 0.5.

Any ideas on this one?


texts: Writing to private dataset in using CARTO's SQL API?

texts:
I know applications can only read from public datasets but I am interested in writing to CARTO from a data collection app to a private dataset.

Is this possible?


texts:
I am making an routing web application similar to Google Maps Get Direction feature. But because of the limited control and data on Google Maps, I decided to use OpenStreetMap instead. This also made me install PostGIS-PostgreSQl database, added routing functionality using pgRouting and now, I'm planning to use a GeoServer. All of these are done in Ubuntu Linux 9.04. My problem is I have difficulty using Linux and I'm more comfortable with Windows. As much as possible, I want to use Windows, but some of the things that I need are tested only in Linux (like osm2pgrouting, which I need to import OSM data to my database).

Is there a way I can implement the project using Windows, to avoid using Linux?


texts:
Which is the best source of free map icons (royalty free & price free)?

I need an icon set for my Google Maps project that will give it a custom feel. Ideally I want an icon set with a creative commons license so I can modify them myself.


texts: What technique would I use to translate any GPS coordinate to a grid?

texts:
If I wanted to divide the existing GPS coordinate system (lat,long) into a "grid" - effectively reducing its resolution - are there any established techniques, formulas or gotchas that I might benefit from knowing about to ensure I'm always translating coordinates to the same grid?

Example:

0.0000000,0.0000000 maps to cell 0,0
0.0686645,0.8047485 maps to cell 1,0

My assumption is that my only input would be a desired size for each cell and likely the assumption that the grid will always originate from the same point (as my example above implies, let's say 0.000000,0.000000).

Moreover, is there a way to do this without having to perform lookups against a large table of cells mapped to GPS coordinate ranges?  Ideally I'd like to be able to just take a latitude and longitude, run it through a formula and know what cell it is within.

Apologies for the very beginner nature of this question.


texts: how do I position a popup relative to a point in Openlayers?

texts:
I have a fixed size FramedCloud Popup that I want to position relative to the point calling it, always above the point and with its right side 50px to the right of the point.

I am looking at the API and see positionBlocks and blocks, but I do not understand how to use them and if that's what I need. I'm really looking for.

I tried setting popup.relativePosition, hoping to find an offset parameter somewhere, but relative position doesn't seem to be having any effect.

var vector = new OpenLayers.Layer.Vector("Points",{
    eventListeners:{
        'featureselected':function(evt){
            var feature = evt.feature;
            var popup = new OpenLayers.Popup.FramedCloud("popup",
                OpenLayers.LonLat.fromString(feature.geometry.toShortString()),
                new OpenLayers.Size(275,71),
                "

" + feature.attributes.title +"

",
                null,
                true
            );
            popup.imageSrc = 'img/popup.png';
            popup.autoSize = false;
            popup.imageSize = new OpenLayers.Size(275,71);
            popup.relativePosition = "bl"
            popup.fixedRelativePosition = true;
            //popup.positionBlocks = {}
            feature.popup = popup;
            map.addPopup(popup);
        },
        'featureunselected':function(evt){
            var feature = evt.feature;
            map.removePopup(feature.popup);
            feature.popup.destroy();
            feature.popup = null;
        }
    }
});


texts: QGIS- how to prevent layers from appearing in legend on composer

texts: My laptop (Lenovo x1) usually crashes When I add legend to a map composer in QGIS (QGIS 3.2 BONN). Problem is that all layers in the project will be "auto-updated" to the legend and this will get my laptop "stuck". Why is this "auto-update" box is pre-defined and checked. Of course this is dependant on your project size..but is there anyway to prevent this? i.e How to code "un-checked" auto-update box from the start?


texts: Association shapefiles to custom data in R

texts:
I have an Australia shapefile (can be found here http://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/1259.0.30.001July%202011?OpenDocument)  and some custom data which contains latitudes and longitudes for a number of people.

How can I join the shapefile to the custom data using latitude and longitude so that I can create a population distribution map and colour the segments of the map accordingly?


texts: QGIS 3.4.11: Labels disapear from map when exporting as picture from the print composer

texts: My question is basicly how to solve the issue mentioned in the topic. I have a map, where labels (city names) are visible both in the project and in the print composer. But when I export the map as a picture, the labels disapear. Is this a bug or do I need to choose some option?


texts: Extracting only elevation data from NASA's GLAS Satellite LiDAR to create raster DEM?

texts:
I'm trying to open data downloaded from NASA's Earthdata site in ArcGIS 10.4.1. Specifically I'm trying to open "GLAS/ICESat L2 Global Land Surface Altimetry Data (HDF5) V034" (Image below, highlighted in blue, can be found by searching 'Iran Glas').

I have the data downloaded, there are multiple .H5 files which I've tried adding to an open ArcGIS .MXD by using the 'Add Data' tool. Once the files are open in the 'Add Data' wizard, there are even more files to chose from, most of these I'm not interested in. It is difficult/near impossible for me to figure out which files are actually important for importing.

I'd like to know if there is an easy way to filter the data so I can extract only elevation data and then, produce a DEM from the extracted data. I've tried opening several files that prompt to produce pyramids for them - this gives me some hope that the data I want is in specific files. However, it is taking an extremely long time for them to open and build on my machine (too long considering how many files there are to open).

I've started looking at converting the .H5 files to .LAS files, mainly prompted by following ESRI's advice (here: http://www.esri.com/esri-news/arcuser/summer-2013/5-ways-to-use-lidar-more-efficiently).

Would this significantly reduce processing time?

I have a reasonable idea of what I want but need advice on the process and whether this is feasible.


texts: How can I bring a layer to the back in QGIS?

texts:
I want to bring to the back (wrap behind) a layer in QGIS and I have no idea how. Anybody knows how?

the layers I want to bring to the front

if I select every layer, the big one pops over


texts: How to organize random spatial data for an area of interest?

texts:
I'm trying to take a huge amount of spatial data for a specific 200,000 acre area, collected over 15 years dating back to 2000, and organize it. The data layers relate to wildfires, habitat rehabilitation efforts, biology, vegetation etc.

The goal is that for any point (or tract of land) within this area, the vegetation/monitoring/wildfire/etc history could be easily accessed, without having a unwieldy mxd full of layers to navigate through.

Is there a method or general framework to organize large amounts of spatial data for a particular area of interest?

I am using ArcGIS 10.2.2


texts: Save QGIS project to PostgreSQL/PostGIS database using PyQGIS

texts:
I'm currently working on some code to automate the creation of a QGIS project from tables in a PostgreSQL DB. For info I'm working in QGIS 3.4.

I'm almost there, but I'm stuck on one thing

How do I save the project back into postgreSQL?

The code I've put together so far:

# CONNEXION
# This works fine

    uri = QgsDataSourceUri()
    DB = "MCD_1"
    user = "user1"
    MP = "XXXX"
    host = "XXX.XXX.XX.XXX"
    port = "XXXX"
    schema = "XXXX"

    uri.setConnection(host, port, DB, user, MP)

# OUTPUT FOLDER
# Here I define the ouput options - saving to local works ok but I'm not sure how to modify the 'PG_OUTPUT' to save into the DB

    LOC_OUTPUT = 'C:/Users/user1/Docs/SIG/test_output/loc_ouput.qgz'
    PG_OUTPUT = 'uri.setConnection(host, port, DB, user, MP)/TEST.qgz'

# LOAD LAYERS
# works ok

    uri.setDataSource(schema, "DB_table_A", "geom")
    layer = QgsVectorLayer(uri.uri(False), "SIG_table_A", "postgres")
    QgsProject.instance().addMapLayer(layer)

    uri.setDataSource(schema, "DB_view_A", "geom","", cb_PrimaryKey)
    layer = QgsVectorLayer(uri.uri(False), "SIG_view_A", "postgres")
    QgsProject.instance().addMapLayer(layer)

# CREATE AND PLACE LAYERS IN GROUPS
# works ok

    root = QgsProject.instance().layerTreeRoot()

    Group_1="Tables"
    group = root.addGroup(Group_1)

    table_1 = QgsProject.instance().mapLayersByName("SIG_table_A")[0]
    root = QgsProject.instance().layerTreeRoot()
    layer = root.findLayer(table_1.id())
    clone = layer.clone()
    group.insertChildNode(0, clone)
    root.removeChildNode(layer)

    Group_2="Views"
    group = root.addGroup(Group_2)

    view_1 = QgsProject.instance().mapLayersByName("SIG_view_A")[0]
    root = QgsProject.instance().layerTreeRoot()
    layer = root.findLayer(view_1.id())
    clone = layer.clone()
    group.insertChildNode(0, clone)
    root.removeChildNode(layer)

## SAVING PROJECT
# Here I just modify the variable in project.write() between PG_OUTPUT and  LOC_OUTPUT - to save locally or in the DB.

    project = QgsProject.instance()
    project.write(PG_OUTPUT)

sources:

https://docs.qgis.org/3.4/pdf/en/QGIS-3.4-PyQGISDeveloperCookbook-en.pdf

Creating Proj and layers:

QGIS Python script for creating project file

https://gis.stackexchange.com/search?q=grouping+layers+pyqgis

Grouping layers using PyQGIS makes them not visible on canvas?

Adding layer to group in layers panel using PyQGIS?


texts: Is it possible to interrupt Select Layer By Location when first feature selected?

texts:
I´ve got a point feature class with some 7.000.000 points and I´m trying to check if any of these points is simply intersected by several polygon feature classes using 'for' loop and 'arcpy.SelectLayerByLocation'. Polygon feature classes contain just 1 polygon feature each and are saved in a folder as shapefiles.

I can create a feature layer from points and then use Select Layer by Location with "INTERSECT" parameter but my problem is that the polygon feature can be intersected by e.g. 700.000 points and the process of selection takes ages.

My question is: is there any way to break the Select Layer by Location process after it finds a first point that is intersected by the polygon? Or do I have to wait until the whole process finishes? The code is simple:

# Create a layer from feature class with 7.000.000 points
arcpy.MakeFeatureLayer_management(pointFC, "point_lyr")

# Builds a list of polygon feature classes within a workspace
fcs = arcpy.ListFeatureClasses("", "polygon")
for polygon in fcs:
    arcpy.SelectLayerByLocation_management("point_lyr", "INTERSECT", polygon, "", "NEW_SELECTION")

    # And now instead of:

    # Count selected points after the process ends and do something if condition is met
    points_count = int(arcpy.GetCount_management("point_lyr").getOutput(0))
    if points_count != 0:
        something

    # I´d rather have

    if a point intersects the polygon for the first time:
        tell me that this happened, then break this loop and start another polygon loop

Is it possible to break a running arcpy process?

Running on Win7, ArcGIS for Desktop Basic 10.1 + Python 2.7


texts: I have contour lines with height values from 0 to 1030 (sorted by height 110, 115, 120, 125). Values with in-between number (like 5, 15, 115, 135...) should be grouped and colored with same color, and other values should be grouped this way 10, 60, 110, 160, 210, 260... (yellow), 20, 70, 120, 170, 220... (red), 30, 80, 130, 180, 230... (cyan) etc. I hope you get it despite my bad English :) So, I don't know how to change symbology in ArcMap to color them that way. Can somebody help me?


texts: Compiling custom plugins for QGIS on Windows

texts:
I've made my first plugin using plugin builder and I want to compile it. I read that I need two commands to do it under windows: pyuic4 and pyrcc4. But when I use any of these commands I get this:

D:\Programy\QGis\bin>pyuic4 -o C:\Users\noaal\.qgis2\python\plugins\Test\Test.py
C:\Users\noaal\.qgis2\python\plugins\Test\Test_dialog_base.ui
ImportError: No module named site

I also read that I must change system path in windows, but I don't know to which location. So anyone have any clue?


texts: Cannot set input into parameter error on second iteration

texts:
On the second iteration of the loop starting line 7 I get the following error.  Everything works as intended on the first iteration.  Any ideas on why this is failing the second time around?

Error:

  Runtime error  Traceback (most recent call last):   File "",
  line 10, in    File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 7221, in
  SelectLayerByAttribute     raise e ExecuteError: ERROR 000622: Failed
  to execute (Select Layer By Attribute). Parameters are not valid.
  ERROR 000628: Cannot set input into parameter in_layer_or_view.

Code:

mxd = arcpy.mapping.MapDocument("CURRENT")
ObjIDs = []
with arcpy.da.SearchCursor('sLineGravityMain_L', 'OBJECTID') as SCur:
    for SRow in SCur:
        ObjIDs.append(SRow[0])

for ObjID in ObjIDs:
    tmpLyr = 'temp' + str(ObjID)
    arcpy.MakeFeatureLayer_management(PathToFC,tmpLyr,'"OBJECTID"=%d' % ObjID)
    arcpy.SelectLayerByAttribute_management('sLineGravityMain_L',"NEW_SELECTION", '"OBJECTID"=%d' % ObjID)
    arcpy.SelectLayerByLocation_management('sManhole_P', "INTERSECT",tmpLyr,"","NEW_SELECTION","")
    arcpy.Delete_management(tmpLyr, 'GPFeatureLayer')

FYI, PathToFC is a placeholder, I have the full path passed in for the actual code. The in_layer for MakeFeatureLayer_management is the 'sLineGravityMain_L' feature class.


texts: GeoServer+GeoWebCache: request does not align to grid(s) 'EPSG:3407'

texts:
I'm having a problem integrating GeoServer and GeoWebCache. When requesting the tiles, the geocaching result is always unsuccessful, and the error message I get is "geowebcache-miss-reason: request does not align to grid(s) 'EPSG:3407'.
I configured the gridset as in this

and the layer bounds
.
If I generate the tiles manually using the Seeding tool, the generated tiles in the preview are always empty. If I use the preconfigured gridsets, it works as expected.
Any hints what I might be missing?


texts: Why does ArcGIS Resample tool change cell size and numbers of columns and rows?

texts: I am working with images of oxygen distribution within individual plant root systems submerged in water. The program I use to acquire these images, LookatRGB, splits the green band in half and makes two green rasters. However, in this process, it changes the cell size to (1.4580744, 325.25). When I use the resample tool in ArcMap, the rows and columns change from the correct number of 1953, 1301 to an incorrect number of 1953, 290212. I have tried changing everything in the environment setting to the correct numbers and projecting to a UTM coordinate system, with no results. As of now, the only way I have been able to fix the problem is by changing the properties in Imagej, then exporting as a tiff. I used this corrected layer in the environment settings, but I need to be able to process thousands of images in the future using ArcMap's model builder/python. Fixing individually in ImageJ then importing into ArcMap will not work.


texts:
I want to upload a csv file with co-ordinates. The locations belong to USA. So I Upload a csv file to ArcMap and add XY data. The Projected Coordinate System is WGS_1984_World_Mercator. As soon as I add a base map, the locations are shown below continent of Africa. I tried changing the co-ordinate system. But not lucky ever. I am a beginner.

Anything I am overlooking?


texts: Controling or modifying FID order using arcpy.AddXY_management

texts:
Using arcpy with slightly modified code from created by Felix Creating parallel line in ArcGIS Pro using ArcPy I made four parallel lines. I then used arcpy.AddXY_management to extract the vertices from the four lines. Up to here, all is good. My problem is controlling or modifying the sequence of which arcpy.AddXY_management assigns FID. Below is a simple version.
Here is a subset of the full version. .

It is hard to depict and explain, but I want to be able to control which point is FID=0 and have sequence proceeded snake-wise. For example, start at the upper right point go down the line to the next point, move one point left, go up to the point across the line, move left to the next point down the line to the next point, and then left again ( lower left point) and then back up across the line. The goal is to create feature points with an FID order for like a lawnmower or a tractor.


texts: "Update Altitude" does not work on TCX Converter

texts:
I have downloaded TCX Converter from the official site.

Although I update altitude data of points which I created from Google Earth, they can't be displayed on table. Only "ml>" is written on altitude(ALT) field instead of integer values.

Why does it happen?(Problematic area is shown in red frame.) ()


texts:
I have a Python application and it uses PyQGIS to - read NetCDF file, convert it to GeoTIFF by writing the rasterlayer with rasterpipe, QgsRasterFileWriter and rasterShader related methods.
However, there is too much pixelization in the TIFF when I zoom in.
input_nc_filename = "input.nc"
rlayer = QgsRasterLayer(input_nc_filename, "east")

resampleFilter = rlayer.resampleFilter()
resampleFilter.setZoomedInResampler(QgsCubicRasterResampler())
resampleFilter.setZoomedOutResampler(QgsCubicRasterResampler())
resampleFilter.setZoomedOutResampler(None)

extent = rlayer.extent()
provider = rlayer.dataProvider()
has_stats = provider.hasStatistics(1, QgsRasterBandStats.All)
stats = provider.bandStatistics(1, QgsRasterBandStats.All, extent, 0)

if has_stats is not False:
    logger.debug("minimumValue = {}".format(stats.minimumValue))
    logger.debug("maximumValue = {}".format(stats.maximumValue))

if stats.minimumValue < 0:
    stats_min = 0
else:
    stats_min = stats.minimumValue
stats_max = stats.maximumValue

width, height = rlayer.width(), rlayer.height()
renderer = rlayer.renderer()
provider = rlayer.dataProvider()

pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())

output_filename = "qgis_cubic.tif"
file_writer = QgsRasterFileWriter(output_filename)

file_writer.writeRaster(pipe, width, height, extent, rlayer.crs())

When I tried setting the re-sampling algorithm to Cubic in QGIS desktop application, it worked perfectly.
Now I am trying to do the same in my Python application. I came across few Python classes related to re-sampling -

QgsRasterResampleFilter
QgsRasterResampler, QgsRasterResamplerV2
QgsBilinearRasterResampler, QgsCubicRasterResampler

Unfortunately, I could not get a breakthrough to use these classes while generating the GeoTIFF.
We use GeoTIFF in a webapp developed in Mapbox GL JS and Angular and not in QGIS hosted environment.


texts: Function to intersect another polygon using OpenLayers

texts:
I wonder if someone here have a JavaScript function that when I draw a polygon I can see if I am drawing on another polygon?

like intersect... if polygon intersect say to me a alert...

using OpenLayers.


texts: GDAL Python Not Working Correctly - Raster Calculations

texts:
I am sure this question has been asked and answered before. I am working with a x64 bit installation of OSGeo4W. I am trying to perform raster calculations using the gdal through the OSGeo4W shell, however every time I try to do the calculation I get the following error:

    C:\>gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B"
Traceback (most recent call last):
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site.py", line 548, in <module>
    main()
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site.py", line 530, in main
    known_paths = addusersitepackages(known_paths)
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site.py", line 266, in addusersitepackages
    user_site = getusersitepackages()
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site.py", line 241, in getusersitepackages
    user_base = getuserbase() # this will also set USER_BASE
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site.py", line 231, in getuserbase
    USER_BASE = get_config_var('userbase')
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\sysconfig.py", line 516, in get_config_var
    return get_config_vars().get(name)
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\sysconfig.py", line 449, in get_config_vars
    import re
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\re.py", line 105, in <module>
    import sre_compile
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\sre_compile.py", line 14, in <module>
    import sre_parse
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\sre_parse.py", line 17, in <module>
    from sre_constants import *
  File "C:\PROGRA~1\QGISVA~1\apps\Python27\lib\sre_constants.py", line 18, in <module>
    from _sre import MAXREPEAT
ImportError: cannot import name MAXREPEAT

I have checked the location for all the required materials, and everything seems to be correct. I have also tried reinstalling the proper gdal python bindings and environmental variables to no avail.


texts: I need shapefiles of Indian states that of 2011 which has almost 671 districts of India. Can anyone please help me?


texts: QGIS Python AttributeError object has no attribute 'pushVersionInfo'

texts:
I am unable to use almost all of the processing toolbox features in QGIS. I keep getting the following Python error across multiple algorithms:

  AlgorithmDialog.py", line 208, in runAlgorithm
  self.feedback.pushVersionInfo() AttributeError:
  'QgsProcessingFeedback' object has no attribute 'pushVersionInfo'

I'm new to QGIS and I am running 3.6.1-Noosa, although I have uninstalled it and updated to 3.6.2 twice. I'm not sure what that is about either.

QGIS code revision 2468226bc9.


texts: Filling NoData gaps in raster using ArcGIS Desktop?

texts:
I have a raster with gaps in it. I want to fill these in with averages of surrounding cells using ArcGIS 10.2 for Desktop.

I have tried "Focal Statistics" but it averages EVERY cell and not just the gaps.

I have tried [Mosaic].IsNull.Con([Average], [Mosaic]) found here Patching
but don't know how to implement it/it won't work for me

I have tried Map ALgebra but I think they were using versions before 10.2.

What I want. Make a fake raster with gaps filled in. Use the fake raster to fill in the gaps of the real raster without replacing the points with actual data.


texts: How to perform IDW interpolation on ocean data points without interpolating on land?

texts:
I am using QGIS 2.14 Essen. I am working on ocean data points and I would like to interpolate them to create a raster. The problem is, when I do the interpolation, the areas where there is land are also filled with data. What I want to do is to interpolate the data only on the oceans area and not on land. I want the land areas to remain with no data.

Is there anyway to this?


texts:
I want to create some transects and a grid with specific lengths, but I realized Qgis is under measuring distances. For the transects I used "Advanced digitizing" (100m lines) and for the grid I used "create grid" with horizontal and vertical spacing of 100 m. However, instead of 100m I got in both cases 75m.
All layers are with the same projection > EPSG:3857 - WGS 84 / Pseudo-Mercator.

Any ideas why this is happening?


texts: Converting a PDF to an Elevation GeoTIFF with QGIS to be imported into Drone Deploy

texts: Is it possible to convert a PDF drawing of a final grade plan, using QGIS into an elevation GeoTIFF to be imported into Drone Deploy?


texts:
I  am using GeoServer 2.15 to provide the INSPIRE view service.  I am testing it with the INSPIRE reference validator and the test fails for the GetMap operation. The validator expects the service to throw an exception when the parameters, “version”, styles”, “crs” and “transparent” are missing. Unfortunately GeoServer does not throw the exception. Looking at the standard specification, it seems that ¾ of these parameters, except “transparent”, are mandatory. I tried to click the “Strict CITE compliance” check box in the GeoServer WMS settings page, but with no luck.

Are there any suggestions on how I could solve this?


texts:
I am adding spatial objects to a GeoPackage database. To register a new spatial table I use a formulation like:

INSERT INTO gpkg_contents
        (table_name, data_type, identifier, description, last_change,
         min_x, max_x, min_y, max_y, srs_id)
VALUES ('location',
        'features',
        'Locations',
        'Relevant locations',
        DateTime('now'),
        -180, 180, -90, 90,
        4326);

I use the validate_gpkg.py utility from GDAL to validate the resulting GeoPackage, and always get the same error:

Traceback (most recent call last):
  File "validate_gpkg.py", line 1673, in <module>
    verbose=verbose)
  File "validate_gpkg.py", line 1640, in check
    checker.check()
  File "validate_gpkg.py", line 1614, in check
    self._check_gpkg_contents(c)
  File "validate_gpkg.py", line 312, in _check_gpkg_contents
    (last_change, table_name))
  File "validate_gpkg.py", line 88, in _assert
    raise GPKGCheckException('Req %s: %s' % (str(req), msg))
__main__.GPKGCheckException: Req 15: last_change = 2019-03-27 16:20:27 for table_name = profile is invalid datetime

Instead of DateTime('now') I have also tried with CURRENT_TIME and CURRENT_TIMESTAMP, but the result is always the same. How can I obtain a GeoPackage-valid time stamp to insert in the gpkg_contents table?


texts: Carto SQL fails due to CORS error in local CartoDB?

texts:
I have successfully installed CartoDB in Ubuntu 14.04 LTS and uploaded a point geometry dataset. I have created a HTML/JS application to access the layers from the local cartodb and visualize in the leaflet client application using cartodb.createLayer.

I need to display a chart by using the data from the uploaded layer. When I try using CartoDB SQL API, I get the CORS error. The details are below.

Can you let me know whether we need to include any configuration for CORS in local cartoDB?

Code:

$.getJSON('http://development.localhost.lan:3000/api/v2/sql?q=SELECT * FROM tablename & api_key = 1a77a491a96dff1c380fc869498587db105c5b65', function(data) {
  console.log(data);
  $.each(data.rows, function(key, val) {
    // do something!
  });
});

Error:

  Cross-Origin Request Blocked: The Same Origin Policy disallows reading
  the remote resource at
  http://development.localhost.lan:3000/api/v2/sql?q=SELECT%20*%20FROM%20%22qgis_view_standard_all_permits_processed_current%22%20&%20api_key%20=%201a77a491a96dff1c380fc869498587db105c5b65.
  (Reason: CORS header ‘Access-Control-Allow-Origin’ missing).


texts: I have a large point cloud (LAS file) and I am trying to thin this using a complex technique. I know it is possible to thin a point cloud using simple techniques such as extracting every nth point. However the point cloud I am working with is more dense in some areas and less dense in other areas. I would like to be able to thin the cloud so that the denser areas are thinned more heavily than the less dense areas. To do this manually one could create a grid and select certain number of points within each cell. (ie one at each corner and one in the center) However this point cloud is much too large for this task. I am trying to figure out if there are any existing solutions to this problem such as a complex selection query, other automated process, or existing software application that could be used.


texts:
I downliaded a SD file from here: http://portal.emodnet-hydrography.eu/geographic-grid
How can I open it in QGIS?


texts:
I have written a python script that iterates through the field of a table using the SearchCursor and produces two table field and the i want to add a number to the individual values in a field produced. The codes seems to work on only one number.Please help what I am missing in my codes:
The code is :

import arcpy
fc = arcpy.GetParameterAsText(2)

rows = arcpy.da.SearchCursor(fc,['S10','S100'])
for i in range(2,7,2):
        arcpy.AddMessage("When i is {0}".format(i))
        for row in rows:
                if row[0]:
                        new_row = row[0]+i
                        arcpy.AddMessage(new_row)

The out put i get:

 When i is 2
   2.275
   2.342
   2.257
   2.667
 When i is 4
 when i is 6

QN: Why is not iterating when i is 4 and 6?


texts: Using ModelBuilder for closest facility analysis in ArcGIS Network Analyst?

texts: I am building a model to iterate closest facility analysis for every building at an archaeological site, with the goal of generating a from-everywhere, to-everywhere aggregate traffic map. I have run the analysis manually on one "incident", and it works perfectly. Now I am trying to build a model to iterate through every building at the site. However, I cannot get the model to generate anything. I press "run entire model" and the dialogue says that it worked, but I can't find any outputs anywhere. I have specified a workspace (folder) and scratch space (folder). Below is a screenshot of the model as currently configured:


texts: What is the difference between != NULL and IS NOT NULL in QGIS filter expressions?

texts:
I have a shape layer with an attribute containing NULL values to which I would like to apply a filter for values different from NULL.

Using the available gui for such query construction, one would intuitively try

"obj_art" != NULL

Which means 'give me all the features with attribute "obj_art" different from NULL' (these count definitely more than 0). Testing this query delivers a strange result from my point of view:

So what I've learned so far is that I can achieve this by using

"obj_art" IS NOT NULL

The question is, what is the difference between != NULL and IS NOT NULL?


texts:
I have a very large dataset with line features that represent a utility, and lines that come off of those utility lines showing connections to buildings.
(example shown below - utility line in black and offshoots in purple.)

I developed some iterative code that takes my selected black lines, selects the offshoots attached to it, and checks some of the attributes. If the attribute checks fail, it takes the id code, appends it to a list, then later uses that list to select the purple lines and export them to their own feature. The code works fine, this is not a troubleshooting post. My questions is about efficiency. Has anyone tried something like this in the past, and is there a better way to do it? I ask because this code took about 10 minutes to run on a given study area, and we have >100,000 study areas. If you multiply that out, it would take ~19 years to run on the whole data set, which essentially renders this useless. The study area had ~130 black lines, and ~450 purple lines in it. The tool took 604 seconds to run, which is a little over a second per feature.

My thought is that running nested for loops that each call upon a Arcpy SearchCursor is very slow, and I was wondering if there way any way to get a similar effect with better performance. Are there any links to documentation or other posts? I assume this is a rather niche use case.

with arcpy.da.SearchCursor("Feature1", ['SHAPE@','ID','Attribute1','Attribute2']) as cursor1:
for row1 in cursor1:
    f_shape = row1[0]
    ID1 = row1[1]
    a1 = row1[2]
    a2 = row1[3]
    arcpy.SelectLayerByLocation_management("Feature2","BOUNDARY_TOUCHES",f_shape)
    with arcpy.da.SearchCursor("Feature2", ['SUBTYPE','ID_1','Attribute3','Attribute4','ID_2']) as cursor2:
        for row2 in cursor2:
            sub_type = row2[0]
            if sub_type <> 1:
                if a1 == 'criteria':
                    if (a2 == 'something1' or a2 == 'something2'):
                        if row2[2] <> 1:
                            found_list.append(row2[4])
                    elif a2 == 'something 3':
                        if row2[3] <> 1:
                            found_list.append(row2[4])
                if row2[1] <> ID1:
                    if row2[4] not in found_list:
                        found_list.append(row2[4])


texts: Seeing all existing joins in any MXD layer in one place?

texts:
I've got an MXD with about fifty layers, and I know there is a table that is joined to at least one, possibly more, of those layers. But I don't know which one :)

I know I can easily look at a specific layer to find what (if any) joins it has. However, is there a way to get an MXD-wide overview of joins that can list all layers with joins, and what they are joined to? (I'm always a fan of shortcuts...)


texts: Delimiting service area for electric network using ArcGIS for Desktop?

texts: I'm working on a study whose objective delimitation of the service area of an electric network with a type polyline I will like to know if there is a way whatever its nature, which will allow me to achieve the result below using ArcGIS 10.2 for Desktop?


texts:
I'm currently trying to create a series of subplots within a larger plot, but outside of a nesting smaller plot, and am having trouble with it. This is the code I'm working from. I've also added a picture of what I'm trying to get at the bottom.

Learning QGIS=group
input_layer=vector
size=number 15
subplotsize=number 2
subplots=number 4
squares=output vector
minisquares=output vector
from qgis.core import *
from processing.tools.vector import VectorWriter
import random
# get the input layer and its fields
my_layer = processing.getObject(input_layer)
fields = my_layer.dataProvider().fields()
#get the input layer and its fields
my_layer2 = processing.getObject(input_layer)
fields2 = my_layer2.dataProvider().fields()
# create the output vector writer with the same fields
writer = VectorWriter(squares, None, fields, QGis.WKBPolygon, my_layer.crs())
writer2 = VectorWriter(minisquares, None, fields2, QGis. WKBPolygon, my_layer2.crs())
# create output features
feat = QgsFeature()
feat2 = QgsFeature()
sx=subplots
sps=subplotsize/2
for input_feature in my_layer.getFeatures():
    # copy attributes from the input point feature
    attributes = input_feature.attributes()
    feat.setAttributes(attributes)
    # create square polygons
    point = input_feature.geometry().asPoint()
    xmin = point.x() - size/2
    ymin = point.y() - size/2
    for tx in range(0, sx):
        px=random.uniform(1, 10) - 5 + point.x()
        py=random.uniform(1, 10) - 5 + point.y()
        print px, py
        square = QgsRectangle(px-sps, py-sps, px+sps, py+sps)
        feat2.setGeometry(QgsGeometry.fromRect(square))
        writer2.addFeature(feat2)
    square = QgsRectangle(xmin,ymin,xmin+size,ymin+size)
    feat.setGeometry(QgsGeometry.fromRect(square))
    writer.addFeature(feat)
del writer

I know this is probably a basic question, but I'm new to coding and all I've managed to do is move around the coordinates of the subplots.


texts: Is there a way to apply the "Select features by rectangle\polygon...." function on multiple layers simultaneously?

texts: I have put a map together using openstreemap data that is intended to be used as a base for different projects. However due to the amount of information (number of layers) the loading time is rather slow. I was wondering if there is a way to select items in all layers within an area by using the selected features only once (well can be a few more times, but not once for each layer) and save them separately, so that I can load only the information within a specified suburb in a project, rather than the entire city?


texts: Add Attribute field and change value in PyQGIS

texts:
I'm building my first script purposed to :

Create fields if they doesn't already exist
Fill value in this fields

The code works, but I have to perform it 2 times when the case requierts to Add Fields (when they did not already exist). In fact, when I launch the code the first time, it adds fields missing but it doesn't fill fields and I get an error message :

  The fields -1 of the feature n°... doesn't exist.

Fields are now built and if I re launch the script, I do not get error anymore and it fills the fields.

How can I fix my script to enable it to perform both parts at the first launch?

from PIL import Image
from PIL.ExifTags import TAGS, GPSTAGS
from xml.etree import ElementTree as ET

layer = iface.activeLayer()
dpr = layer.dataProvider()

apr_ind = layer.fields().indexFromName("apr")

# add fields
fields_arr =  [("apr", QVariant.Double)]
flds = [QgsField(i[0], i[1]) for i in fields_arr if i[0] not in [fld.name() for fld in layer.fields()]]
layer.dataProvider().addAttributes(flds)
layer.updateFields()
# field update
layer.updateFields()
layer.startEditing()
for feature in layer.getFeatures():
    file_path = feature["Path"]
    print(file_path)
    img = Image.open(file_path)
    exif = { TAGS[k]: v for k, v in img._getexif().items() if k in TAGS }
    with open(file_path, "rb") as f:
        s = str(f.read())
    start = s.find('<x:xmpmeta')
    end   = s.find('</x:xmpmeta')
    xmp   = s[start:end+12].replace("\\n","\n")
    tree  = ET.XML(xmp)
    apr = exif["FNumber"][0]/exif["FNumber"][1]

    dpr.changeAttributeValues({ feature.id():{apr_ind: apr }})
layer.commitChanges()
print("COMPLETED")


texts: Why does the SAGA raster calculator fail when used in a Processing model?

texts:
I have a problem with the modeler and the saga raster calculator. When I'm using the raster caculator without the modeler it works fine. When I'm using the saga raster calculator in my model I always get the messege: "Oooops! The following output layers could not be open test: C:\Users\miller\AppData\Local\Temp\processing\4a78580e485a40d199a2fa138960042c\RESULTALG2.tif" –

Is this a known problem? Some weeks ago everything worked fine.


texts:
I often use the brilliant SQL Server 2008 Spatial Tools to import shapefiles into SQL Server 2008. This works absolutely fine but it leaves me wondering, how do Microsoft expect end users of SQL Server to get shapefiles into SQL Server? What is the official method for importing this data?

The only references to this I've spotted on GIS.stackexchange are these questions: How to import shapefiles into MS SQL 2008 and then view that data using QGIS? and How to import shape file into SQL Server 2008 R2 using only SQL scripts? , which references spatial tools, GDAL etc. but nothing from Microsoft. Is there something I've missed? Call me naive but I'd almost expect there to be an option in the data source list in SQL Server Import/Export wizard for shapefile!

In summary: I'm not asking how to import shapefiles into SQL Server, I'm specifically asking how Microsoft expect us to import shapefiles into SQL Server.


texts:
I'm using OpenLayers 3.6.0. I want to use ol3-cesium. Ol3-cesium use cesium for glob visualization. I have a question:

In ol3-cesium, Can I see Ups and downs of land? In other words Does it accept DEM files?


texts: Speeding up script to flag nodes within 100 feet of segments using PostGIS?

texts:
I'm working with two datasets, attempting to flag when an intersection is within 100 feet of a network subset. There are over 40,000 intersections and my network is ~ 1400 segments.

The problem with the SQL script below (that I think works), is that the runtime is far too long. I'm fairly new to PostGIS and feel there is likely a faster way to do this. Basically, I'd like the each intersection (row) to have a 1 or 0 where 1 indicates the presence of a network segment within 100 feet.

Anyway to speed this up?

UPDATE intersections g
SET facility_flag =
      CASE
      WHEN EXISTS(
        SELECT 1
        FROM network_subset d
        WHERE ST_INTERSECTS(g.geom,
          ST_BUFFER(d.geom, 100))
      )
        THEN 1
        ELSE 0
      END
    ;


texts: Converting Lat,Long to XY in ArcGIS API for Python?

texts:
I am working on a Alexa skill where I get the user's location (Lat, Long) from the device. Since my skill only applies to a specific city (Boston, USA), I am using the ArcGIS reverse_geocode to determine if that Lat, Long is in Boston, but the input is in the XY format.

I wrote the following method:

from arcgis.gis import GIS
from arcgis.geocoding import reverse_geocode
from arcgis.geometry import from_geo_coordinate_string

gis = GIS()
location = [42.3937647, -71.1449038]
m_location = reverse_geocode(location)
print(m_location)

# Error:
# Cannot perform query. Invalid query parameters.
# Unable to find address for the specified location.

Therefore, I need to convert my lat,long to XY. I found the find_transformation function, but don't know how to use as inputs (the in_sr and out_sr, for instance).

I cannot use packages such as pyproj, because the team does not approve, and I'm positive that ArcGIS has a function for converting between coordinates


texts: batch save as layer to relative path in ArcPy

texts:
I would like to use a similar code as below, but one that uses relative paths. Since my input is a batch of tif files (not layers), I do not know if I need to use other functions like arcpy.SaveToLayerFile_management or os.path.

import arcpy

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

for lyr in arcpy.mapping.ListLayers(mxd):

    lyr.saveACopy(lyr.name + ".lyr")


texts: Discovering highest point in area using QGIS?

texts:
I'll preface by saying I've never used or even heard of GIS software previously, so I'm sorry for being very clueless.

My goal is to find the highest point in Highgate (in London, England).

I have downloaded the "OS Terrain 50" from the Ordinance Survey download page, and have opened it up in QGIS. Unfortunately, I think I misunderstood as this only shows me the contours.

Where should I be looking for height data, and how would I interpret it? I expect I'm struggling because I'm Googling all the wrong things.


texts: Ability to cache geographic data and only update style in Geoserver

texts:
I have a setup where I have a SQL Server database with two tables: one table storing the geometry of the shapes I want to display, another table with some attribute value that is used to color and label the shapes. Then I use this data to generate WMS tiles that are displayed in a web front-end to the end user. In my data, the geometry does not change, but the attributes change very often.

In SQL Server, I configured a view joining the two tables, then used this view as a data source for Geoserver. Then I configured a SLD using these attributes to get the style I want.

It works. However, I have some performance issues, which seems to be linked to the necessity of reloading the geometry from the database for each query. My first idea was tile caching, but I can't use it because there are so many possible values for the attributes, and they change very often.

Ideally, what I would like to do is to have Geoserver cache the geometry data and for each new WMS query only look up the attributes to use and generate a style based on that, instead of retrieving both geometries and attributes at the same time. Is it possible to achieve this?

I am currently using Geoserver, but if it is possible to do it with a different WMS server, I would be interested as well.


texts: Creation of polygon in the gap between two existing polygons

texts: If there is a gap between two polygons is it possible to fill it? Much like what happens using the hatch function in AutoCAD. The output would be a new polygon in place of the gap.


texts: How to use sum of circular moving window for each center cell with focal in R

texts:
I want to calculate the sum of a moving window circle in R for each center cell. Is there a way to do this in R with focal or focalWeights on my raster (preferrably with map units)? Focal does not take circle or map units and focalWeight gives only weighted sums.

I am looking for an alternative to FocalStatisticas (in Esri) in R.

I found a thread at http://r-sig-geo.2731867.n2.nabble.com/Raster-package-Focal-sum-in-circles-td7587683.html which is tackeling quite a similar idea, which made me come up with a code similar to this.

fwModel <- focalWeight(modell_soil, 25, type='circle')
fwModel[fwModel>0] <- 1
p_focal<-focal(modell_soil, w=fwModel ,fun=sum , na.rm=TRUE)

This code allowes me to use the focal command with a circluar window and map units. The only problem I have until now, is that it does not allow me yet to only calculate the sum for the center cell instead of the whole circular window. I tried quite some things but did not get closer to the solution... Is there a way to tackle the center cell that could fit in my method above?

NAs are actually not a problem for me, as I do not have any in my raster. It's just part of the thread (in the link, so I included it in my trie-outs).

When I had to use rectangular windows I already learned to tackle the center cell, but I am not sure how to use it with a circular window.


texts: Use multiple double attributes to symbolize a line shape in ArcGIS

texts:
I need to represent multiple attributes to a line shape.

The attribute table already have the fields I need to represent. One of the attributes is the total number of commuters between cities; this must be symbolized by graduated symbols. The other one is the proportion of commuters in relation to the total population of the city of origin; I would like to symbolize this by graduated colors. The result, specially the legend, should look alike this picture:

Is there a way to do this on ArcGIS?

Ps.: using multiple attributes I figure out that the representation of the quantity by category is not the solution of the question, because what I need is something like "quantity by quantity", moreover, even if the quantity by category was the answer, the problem with the legend would remain unsolved.


texts:
Can you help me with WFS on GeoExt?
I'm read tutorial and see this code.

proxy: new GeoExt.data.ProtocolProxy({
        protocol: new OpenLayers.Protocol.WFS({
            url: "/geoserver/ows",
            version: "1.1.0",
            featureType: "parks",
            featureNS: "http://medford.opengeo.org",
            srsName: "EPSG:4326"
        })
    })

After this I see this var vectorLayer = new OpenLayers.Layer.Vector("Editable features"); And don't understand how to create a layer. So I have some questions.

How to get URL for the WFS? (I use GeoServer)
What does it mean by featureType option.
And more important. This code is correct? I find it here http://workshops.opengeo.org/geoext/wfs/grid.html


texts: Conditionally labeling layer based on three fields

texts:
I am labeling our power line span layer based on the conductor/wire type.

This information is located across three fields in the span layer which represent phases A,B,C.

Basically I want code that will:

Display conductor_a,

if conductor_a is blank display conductor_b,

if conductor_b is blank display conductor_c

I am unsure if python or VB would be better for this. Unfortunately I am having difficulty with entering python into the provided window. I've managed to make the following VB code

[gs_conductor_a] & vbCrLf & [gs_conductor_b] & vbCrLf &  [gs_conductor_c]

This displays the conductor name, but might display up to three labels, which creates a large visual mess in some instances. I want to create a labeling system where there is only one label displayed per line segment. It is worth noting if more than one phase is populated the conductor is the same for all populated fields so anyone of the three will do, but it needs to display only one label total.


texts: How do I obtain the distances of each point lying inside or outside buffers from the center of buffers?

texts:
My aim is to obtain the shortest distance between all the points overlain on a polygon and the center of buffers also overlain on the same polygon.

I have done the following to illustrate this:

import pandas as pd
import geopandas
import numpy as np

df = pd.DataFrame(
       {'City': ['Fargo', 'Orange', 'Jersey City'],
        'State': ['Nevada', 'California', 'New Jersey'],
        'Latitude': [46.877186  , 33.787914 , 40.728157 ],
        'Longitude': [-96.789803, -117.853104, -74.077644]})
df['Coordinates']  = list(zip(df.Longitude, df.Latitude))
df['Coordinates'] = df['Coordinates'].apply(Point)
gdf = geopandas.GeoDataFrame(df, geometry='Coordinates')
gdf_points=gdf
print(gdf.head())

world= geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
#nybb=  geopandas.read_file(geopandas.datasets.get_path('nybb'))

# We restrict to USA.
ax = world[world.name == 'United States'].plot(
    color='white', edgecolor='black',figsize=(20,10))
us_states.plot(ax=ax,figsize=(20,10))

# We can now plot our GeoDataFrame.
gdf.plot(ax=ax, color='red')
ax.set(xlim=(-140,-50),ylim=(0,60))

us_name = world[world.name=='United States']
gdf['Coordinates']=gdf.buffer(3)
#gdf.head(5)
gplot=geopandas.overlay(us_name,gdf,how='identity')
gx=gplot.plot(edgecolor='k', alpha=0.5, cmap='tab10', ax=ax)
gdf.plot(ax=ax)

df1=pd.DataFrame(
      {'Points':['Point1','Point2','Point3','Point4','Point5','Point6'],
         'Latitude':[44.773491,30.287765,41.279906,45.234565,31.23345,38.778987],
         'Longitude':[-94.789323,-114.565309,-72.055234,-93.230988,-112.233456,-72.239989]})
df1['Coordinates']  = list(zip(df1.Longitude, df1.Latitude))
df1['Coordinates'] = df1['Coordinates'].apply(Point)
gdf1 = geopandas.GeoDataFrame(df1, geometry='Coordinates')
print(gdf1.head())

gdf1.plot(ax=gx,color='red')
gx.set(xlim=(-140,-50),ylim=(0,60))
centroid_of_buffers=gdf['Coordinates'].centroid
min_dist=centroid_of_buffers.distance(gdf1[gdf1.Coordinates=="POINT (-93.230988 45.234565)"])#for point4 this is the coordinate.Interpret as the distance between the center of each buffer and point4

But this shows no result for each buffer center distance between it and point4 is NaN


texts: Creating windows form based tools as add-ins in ArcMap using vb2010

texts:
I am writing a project in vb2010 in which i want to add a package having 4 toolbars with about 48 tools in them in ArcMap 10.0.
So far, I have created an ArcMap class project with the common functions in it. I have created windows forms and referenced this function class to it. I ahve also created a ArcMap add-in in which i have created a tolbar and a button.

Now, I want to add this form as a tool in the toolbar on clicking on the button in the add-in.

Is there any other method to achieve the same?
How to add the windows form as a tool in ArcMap Toolbar add-in?


texts: Unwanted Italicization of Symbol Text (ArcMap 9.3.1)?

texts:
So I am running in to a problem that has got me stumped. I'm working on a road systems map, and there are a few symbols with text in the middle. The problem is that the text within the badge, triangle, etc. is randomly becoming italicized. I have gone into the symbol property editor, made sure that the text layer has no italics, rotation, anything. However the text is still being represented in italics when on the map? The symbol itself is not rotated or shifted in any way, just the text inside. I might be missing something simple, but I cannot figure it out for the life of me.

I am using ArcMap 9.3.1, and if it matters, the MPS Atlas Extension.


texts: Clustering (ConvexHull) data based on distance

texts:
First I would like to mention that I am not an expert with PostGreSQL so my question may seem dumb, but here it goes:

I have a database with over 800 000 points associated to a geometry (EPSG 2950) and I would like to cluster those in convex hulls based on this:

"Group every point that has (at least) 3 other points at a distance of maximum 100 meters."

The resulting laye would show polygons delimiting points that respect that condition.

I am guessing I would use the ST_ConvexHull fonction for this particular request but I am struggeling to see how i add my constraints.


texts: using existing SLD style with postgis layers in Opengeo Composer

texts:
I have exported some styles from QGIS as SLD files and then loaded them into GeoServer.  I have published a couple of layers from a PostGIS database and these load nicely in the OpenGeo Composer with a simple YSLD style.

I now want to apply the more complete style that I built in QGIS so I went back to Geoserver and found 'available style' under the 'publish' tab of the layer.  I added the style that I had imported to the 'available' list but composer appears to ignore it.

It seems that the only way I can fiddle with the style is to edit the yaml.  I did try an XML to YAML converter and was unsurprised when the results were not usable.

I really want to have one definitive source for my styles. I don't mind exporting them from QGIS as files but I really don't want to have to recreate them in YAML by hand.

What am I missing?


texts: Creating a really large shapefile without eating all the virtual memory

texts:
I am doing something similar to below to create a really large (large than memory on my machine) shapefile:

shapes = shapefile.Writer(shapefile.POINT)
shapes.autoBalance = 1
shapes.field('MyId','N')
shapes.field('OtherData','C', 20)

for item in largeIteration():
    . . . . .
    shapes.record(item['id'], item['desc'])
    shapes.point(item['lat'], item['lon'])

The problem is while I'm using yield in largeIteration() to not hold onto memory there, the entire shapefile seems to remain open in memory, and therefore shapes grows with each call to shapes.record() and shapes.point(). Is there any way I can have pyshape not hold the entire shape in memory? I see OGR and its virtual file systems. However, I don't see a "stream to disk" vfs option.


texts:
I am interested in exploring the idea of generating a predictive model for a target variable in a spatial context. I would then like to correct the residuals of this model using kriging, i.e. in a nutshell you sum up the results from both models and take that as the final prediction:

The classical approach is to fit a GLM and then krige (is that how you turn it into a verb?) its residuals. I would like to fit something like a Random Forest. Which leads me to my question. How do you validate these models? If I wanted to do cross-validation do I separate the data (e.g. 10-folds) and reapeatedly fit the model, krige the corresponding residuals -> predict on test set, etc. Or do I have to separate the data into 3 sets, one to fit the model, see how it predicts on another data set, krige THOSE errors and then test on the final data set, repeat. What is the usual apporach to do this correctly?

Thank you in advanced.

I'm using R by the way.


texts: How do you access PostGIS geometry data from .NET

texts:
I have a PostgreSQL / PostGIS database with both geometry and non spatial columns
I created a new Console App .Net framework project  and added these nuget packages
EntityFramework
pgsql
EntityFramework6.npgsql
Npgsql.NetTopologySuite
When I create a database first ADO.net data Entity Modal I can see the  non spatial columns but not the geometry see pics
Could someone explain what I need to enable spatial queries?


texts: Mapinfo appending TAB file error - multiple points appearing in the same location

texts: I have two TAB files with identical table structures and when appending one to the other all appended points get placed in one location, however when I append them the other way round it works fine. Both have identical projections and table structures. Any idea what causes this so I can look out for it in the future?


texts: How to show new points in arcgis online collector map?

texts: I have created Arcgis Online Collector Map. Is it possible to have new data points collected to stand out from the old ones? As of now, they all look the same and its hard to tell which ones are new and which ones are old/non existant.


texts: How to create autoincrementing primary key with dynamic writer in FME

texts:
I'm writing 3 shapefiles to an MSSQL Spatial database with FME. I'm using a writer with dynamic settings to simplify future updates of the input files.

In order to create a spatial index with the dynamic writer I need a primary key column. Right now, I'm accomplishing this with a Counter that writes to a manually created "pid"-attribute marked as primary key, but this feels like a hack.

So, how do I create an auto-incrementing primary key on a dynamic database writer in FME?


texts: Using Soil Line to Assess Soil Fertility?

texts:
I am using a Soil Line model to research possible soil fertility in bare soil patches in southern Africa. I am using a function in R called BSL (Bare soil line) to create the soil line from some satellite image: function is at https://rdrr.io/cran/landsat/man/BSL.html. The Soil Line plots NIR and Red reflectance values from a satellite image.

If I am just generating the soil lines to compare the soil lines' slope and  intercept values to estimate soil organic matter, do I need to analyze the vegetation pixel values as well?

I have included the R-generated Soil line graph as an example:

I am new to studying soils so but I have read much of the pertinent literature, such as these ones:

  A.R. Huete, A Soil Adjusted Vegetation Index (1988)

  Baret, F., Jacquemoud, S., & Hanocq, J. F. About the Soil Line
  Concept in Remote (1993)

  Ahmadian et. al. A New Concept of Soil Line Retrieval from Landsat 8
  Images for Estimating Plant Biophysical Parameters (2016).


texts: Using one shapefile to analyze different plots in GEE

texts:
I have one shapefile of agriculture plots that looks like this-

I want to run analysis in google earth engine but for each plot.

Is there any way to run this analysis without saving each plot and upload it as table?

Can google earth engine run my analysis for each plot separately while they are all in one shapefile?

this is my code-

var geometry=table;//this is the featurecollection

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000)
  .copyProperties(image, ['system:time_start']);
}

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2019-06-01', '2019-06-30')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .select('B1','B2','B3','B4','B8','QA60')
                  .filterBounds(geometry)
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

var clippedCol=dataset.map(function(im){
   return im.clip(geometry);
});

//test if clipping the image collection worked
Map.centerObject(geometry,9);
Map.addLayer(clippedCol.median(), rgbVis, 'RGB');

//function to calculate NDVI
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4'])
  .rename('NDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(ndvi);

};

//NDVI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI).select('NDVI');

var colorizedVis = {
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};

//analyze images from image collection collection
var listOfImages = withNDVI.toList(withNDVI.size());
var listOfNumbers =[0,1,2];

var myImage = withNDVI.first()

var collectData = function(feature){
  var myGeom = feature.geometry()
  var myReduction = myImage.reduceRegion({
                  reducer: ee.Reducer.mean(),
                  geometry: myGeom,
                  scale: 30            // Change to match your need
    })
  return(ee.Feature(myGeom,{mean: myReduction}))
}

var myNewFeatures = geometry.map(collectData) ;

for (var i in listOfNumbers) {
  var image = ee.Image(listOfImages.get(listOfNumbers[i]));
  var meanDictionary = image.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry:myNewFeatures.geometry()
  });

 var STDDictionary = image.reduceRegion({
  reducer: ee.Reducer.stdDev(),
  geometry:myNewFeatures.geometry()

});

// print(i,'mean',meanDictionary,'standard deviation',STDDictionary);

var std2 = ee.Number(STDDictionary.get("NDVI")).divide(2);
var mean1 = ee.Number(meanDictionary.get("NDVI"));

var negBorder=mean1.subtract(std2);
var posBorder=mean1.add(std2);

var imageNDVI=image.select('NDVI');
var gtPOS=imageNDVI.gt(posBorder).selfMask().rename('PositiveBorder');
var ltNEG=imageNDVI.lt(negBorder).selfMask().rename('NegativeBorder');
var betMEAN=imageNDVI.gt(negBorder).and(imageNDVI.lt(posBorder)).selfMask().rename('MeanBorder');

var PositiveCOL = {
  palette: [
  '006622'
  ],
};

var NegativeCOL = {
  palette: [
  'ff0000'
  ],
};

var MeanCOL = {
  palette: [
  '00ff00'
  ],
};

Map.addLayer(gtPOS,PositiveCOL,i);
Map.addLayer(ltNEG,NegativeCOL,i);
Map.addLayer(betMEAN,MeanCOL,i);

 }


texts: MMQGIS - Failure conneting to API (WinError 10061)

texts:
After setting up the the Geocode/Geocode CSV with Google/OpenStreetMap I get the following error:

My .csv is set up in UTF-8; QGIS release is 3.14 Madeira.

Any ideas what is causing the abovementioned error?