texts: Getting working Python script to run from toolbox?

texts:
I'm creating a model, using ArcGIS ModelBuilder, and want to use SelectRandombyPercent from Randomly subsetting X% of selected features using ArcPy? to randomly select points from a feature class.  The selected features will be an intermediate data set for use as a source for a Network Analyst OD Cost Matrix layer destinations.  The script works perfectly from the Python shell in ArcMap.

I've added a Script in my toolbox utilizing the script from the shell.  When executed, I'm prompted for parameters, get the tool executing message, and get a Succeeded (green check) result. However it doesn't actually make a selection.

What do I need to modify to have the script (run in ModelBuilder) make a selection that I can use?

It appears that the major problem with moving from Shell to Script/Tool was the function definition line (probably a rookie mistake, but definitely not easily discoverable).

FINAL WORKING CODE (verified in both toolbox and ModelBuilder):

import arcpy
import random

inLayer = arcpy.GetParameterAsText(0)
inPercent = arcpy.GetParameterAsText(1)

arcpy.AddMessage(inPercent+"% random sample")

if float(inPercent) > 100:
    arcpy.AddError ("Error: Percent is greater than 100")
    sys.exit()

if float(inPercent) <= 0:
    arcpy.AddError ("Error: Percent is less than zero")
    sys.exit()

fc = arcpy.Describe (inLayer).catalogPath
featureCount = float (arcpy.GetCount_management (fc).getOutput (0))
count = int (featureCount * float (inPercent) / float (100))

if not count:
    arcpy.SelectLayerByAttribute_management (inLayer, "CLEAR_SELECTION")

oids = [oid for oid, in arcpy.da.SearchCursor (fc, "OID@")]
oidFldName = arcpy.Describe (inLayer).OIDFieldName
delimOidFld = arcpy.AddFieldDelimiters (inLayer, oidFldName)
randOids = random.sample (oids, count)
oidsStr = ", ".join (map (str, randOids))
sql = "{0} IN ({1})".format (delimOidFld, oidsStr)
arcpy.SelectLayerByAttribute_management (inLayer, "", sql)

arcpy.MakeFeatureLayer_management(inLayer,outLayer)

I made five major changes from the source:

Deletion of the DEF line.
Addition of FLOAT for the percent comparison in the IF lines.
Addition of message with the calculated percentage (for quick reference).
Addition of an error message and function stop in the IF lines.
Addition of MakeFeatureLayer parameter and output.


texts: Unable to get ogr2ogr to recognise ESRI file geodatabase

texts:
I'm trying to get ESRI File Geodatabase working on a Ubuntu 18.04 instance.  I have followed the instructions in the readme contained in https://github.com/Esri/file-geodatabase-api/blob/master/FileGDB_API_1.5.1/FileGDB_API_1_5_1-64gcc.zip

I have GDAL version 3.0.4 installed.

Everything completes, but when I get

ogrinfo --formats | grep GDB

I get the following

OpenFileGDB -vector- (rov): ESRI FileGDB

From my investigation this appears to be down the version of gdal I have installed.  I have tried different versions of gdal from the Ubuntu Official repo and also the ubuntu unstable gis one.

Unfortunately company policy prevents me from compiling packages myself

How can I resolve this?

Edit:

I have now tried compiling this from source

Downloaded https://github.com/Esri/file-geodatabase-api/blob/master/FileGDB_API_1.5/FileGDB_API_1_5_64gcc.zip and extracted it to /opt/FileGDB

Then added export PATH="$PATH:/opt/FileGDB_API"

to  my .profile file
then run source .profile

sudo apt-get build-dep gdal
apt-get source gdal
edited the debian/rules file to add --with-fgdb=/opt/FileGDB_API
dpkg-buildpackage -us -uc

It fails with

/usr/bin/ld: skipping incompatible /opt/FileGDB_API/lib/libFileGDBAPI.so when searching for -lFileGDBAPI
/usr/bin/ld: skipping incompatible //usr/local/lib/libFileGDBAPI.so when searching for -lFileGDBAPI
/usr/bin/ld: cannot find -lFileGDBAPI
/usr/bin/ld: skipping incompatible /opt/FileGDB_API/lib/libfgdbunixrtl.so when searching for -lfgdbunixrtl
/usr/bin/ld: skipping incompatible /opt/FileGDB_API/lib/libfgdbunixrtl.a when searching for -lfgdbunixrtl
/usr/bin/ld: skipping incompatible //usr/local/lib/libfgdbunixrtl.so when searching for -lfgdbunixrtl
/usr/bin/ld: skipping incompatible //usr/local/lib/libfgdbunixrtl.a when searching for -lfgdbunixrtl
/usr/bin/ld: cannot find -lfgdbunixrtl
collect2: error: ld returned 1 exit status
GNUmakefile:55: recipe for target 'libgdal.la' failed
make[3]: *** [libgdal.la] Error 1
make[3]: Leaving directory '/home/stuart/gdal/new/gdal-2.2.3+dfsg'
GNUmakefile:64: recipe for target 'check-lib' failed
make[2]: *** [check-lib] Error 2
make[2]: Leaving directory '/home/stuart/gdal/new/gdal-2.2.3+dfsg'
debian/rules:142: recipe for target 'override_dh_auto_build' failed
make[1]: *** [override_dh_auto_build] Error 2
make[1]: Leaving directory '/home/stuart/gdal/new/gdal-2.2.3+dfsg'
debian/rules:77: recipe for target 'build' failed
make: *** [build] Error 2
```


texts: Why is Spatialite showing dates in MM-DD-YYYY format?

texts:
I have a table with a dates column that is formatted in MM-DD-YYYY format, as shown in SQL SELECT statements.  I'm using Qspatilite (QGIS 1.8) and Spatialite GUI.

When I ORDER BY the date column it sorts by MM, then by DD, then finally by YYYY.
In other words, I get rows like this:

  01/04/2004
  01/24/1987
  02/07/1994
  03/17/1982
  03/27/2009
  ...
  12/29/1997

I understand that SQLITE doesn't enforce data types.  Is there an easy way to convert the format to YYYYMMDD (so that I can sort by date) or do I need to manually parse the date myself?


texts:
I have created a polygon shapefile in Erdas Imagine 9.2 from a raster image. I want to export it to .dxf or .dgn for further processing. When I export the .shp to .dxf in Global Mapper, large polygons sharing the border of an image are missing. I tried in Erdas Imagine but the result was in Line String.

Would appreciate any help.


texts: Is it possible to build a .vrt file from multiple files with Rasterio?

texts:
I would like to build a vrt file from multiple dataset. I know the gdalbuildvrt but i don't find how I'm supposed to do the same using rasterio vrt object.
Does anyone have a simple example ?


texts: Merge Raster Layer Specific Order in Graphic Modeler

texts:
I'm using the Graphic Modeler/Processing Modeler to automate a process. The last step in this process is Merging three different Raster Layers. I'm using the "Merge Raster Layers" from Saga 2.0.8. Is if it is possible to specify the order in which these are merged?

I would like it in this order.

-#1 WaterStreams - Top

-#2 Grid with WaterStreams area (Watershed?)- Middle

-#3 Height Model (DEM) - Bottom

So #3 should be at the bottom, then #2 will be added on top of this and then #3 will be added at the top.

Pretty much a DEM with WaterShed burned into it, with WaterStreams burned through the DEM and WaterShed. If this makes sense?

If I do it manually and use the toolbox and Merge Raster Layers it's no problem (First value in grid order), but in the automated process I can't specify the correct order.

Update:
http://i.imgur.com/phWe2ja.png

Update1::Tried editing the .model file, rearranging the order of #1, #2, #3, but no luck, either I'm not doing it right, or it doesn't work that way.

Update2:Tried a few more things using textpad editing, also tried doing two different "Merge Raster Layer" after each other, but no luck. It always puts #3 at the top, though I need it to be at the bottom.


texts: Calculating distance between points using ArcGIS Desktop?

texts:
I want to write Python syntax to calculate the distance from point to other points in the water network and so on for others, I have 21 points.

I need distances have to follow the network

That means I will have 441 records.

I don’t want to use route method it’s not useful.


texts: Summing flow from cities listed in another layer within QGIS?

texts:
I have transport flow represented by arrows within QGIS in one layer. The database is coded with cities ID for the origin and cities ID for the destination of the flow.

I also have in another layer, transporte lines with id of cities crossed by this line (using this answer to Creating list of intersected polygons for each line using QGIS?)

I would like to get the sum of flow between cities crossed by my transport lines

Line 1 : sum of 'flow' case when 'ORIGIN' is like 'CITIES_INTERSECTED'
AND when 'DESTINATION' is like 'CITIES_INTERSECTED'
Line 2: sum of
'flow' case when 'ORIGIN' is like 'CITIES_INTERSECTED' AND when
'DESTINATION' is like 'CITIES_INTERSECTED'

Something like this i suppose


texts: How can I get all streets width of Paris

texts:
I'm trying to query OpenStreetMap via QGIS (for example) to get a list of all streets in Paris and their width.

example :

Boulevard Sebastopol : 30m
Rue Tiquetonne : 11m
etc

Is there such data in OpenStreetMap ? or anywhere else ?


texts: How to get features based on their attributes using Python - special characters?

texts:
I want to get features based on their attributes and it works fine when I use this Python script:

canvas = iface.mapCanvas()
curLayer = canvas.currentLayer()
it = curLayer.getFeatures(QgsFeatureRequest().setFilterExpression('"ELAVSNAVN" = \'Fraugde By, Fraugde\' and "MATRNR" = \'15g\''))

for baseFeature in it:
    print baseFeature.id()

But nothing is returned when I use special danish letters as shown in this example:

it = curLayer.getFeatures(QgsFeatureRequest().setFilterExpression('"ELAVSNAVN" = \'Åsum By, Åsum\' and "MATRNR" = \'10a\''))

What can I do solve this problem ?


texts: R- Import feature class from geodatabase and plot specific variables

texts:
I have a a big feature class in geodatabase (~1.3 million rows and 234 columns; you could download it using this link) in ArcGIS 10.4.1.

I want to import the feature class from the geodatabase to R and plot specific variables (using ggplot2).

I tried importing it into R, using the answer of this question, as below

require(rgdal)
# The input file geodatabase
fgdb = "C:/Question_Online/union_question.gdb"
# List all feature classes in a file geodatabase
subset(ogrDrivers(), grepl("GDB", name))
fc_list = ogrListLayers(fgdb)
print(fc_list)
# Read the feature class
fc = readOGR(dsn=fgdb,layer="union_thirdtrial")

However, it took more than an hour, so I had to close RStudio to terminate the process.

Can anyone point out what I am missing for a better (faster) way to import this feature class to R?


texts:
Loading *.adf files into QGIS
How do I load .adf files into QGIS?


texts: Error while exporting from Google Earth Engine

texts:
I run an unsupervised classification with Google Earth Engine using this code:
    // UNSUPERVISED classification of Conifers//

var input = ee.ImageCollection("COPERNICUS/S2_SR");

var geometry =
    /* color: #98ff00 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[7.582719257752308, 46.02568420992065],
          [7.582719257752308, 44.97372238809936],
          [9.228269435087903, 44.97372238809936],
          [9.228269435087903, 46.02568420992065]]], null, false);

input = input.filterBounds(geometry);
print(input);

input = input.filterDate('2020-08-06','2020-08-08');
print(input);

input = input.sort('CLOUD_COVER');

input = input.mosaic();
print(input);

input = input.clip(geometry);
print(input);

Map.centerObject(geometry,9);

Map.addLayer(input,{min:-157.33,max:1421.15,bands:['B4','B3','B2']}, 'input image');

var training = input.sample({
  region: geometry,
  scale: 10,
  numPixels: 5000
});

var clusterer = ee.Clusterer.wekaKMeans(15).train(training);

var result =  input.cluster(clusterer);

Map.addLayer(result.randomVisualizer(), {}, 'Unsupervised Classified Image');

Export.image.toDrive({
  image: result,
  description: 'Unsupervised Classified Image',
  scale: 10,
  region: geometry
});

When I run the code with 10m scale in the training part and choose 10m also while exporting, it never finished the export to drive and I got this error:
Error: Export too large: specified 214533809 pixels (max: 100000000). Specify higher maxPixels value if you intend to export a large area.
However, when I set the scale in code and while exporting to 20m, it got exported in 8 minutes.
I need the result in 10m resolution so how can I change the code so it exports a 10m resolution? This is my first attempt in GEE and I am not familiar with how the scale and numPixel have to change.
Please let me know if I can modify any other part of the code.


texts:
I use a python module I created that does tests for point in polygon.
It uses the python osgei/ogr and osr modules.

Specifically, I use the following elements of these modules:

osr.SpatialReference() (with ImportFromEPSG4326 and InputFromEPSG3857)
osr.CoordinateTransformation()
ogr.Geometry() (with ogr.wkbLinearRing, ogr.wkbPoint and ogr.wkbPolygon)

The software using these osr/ogr capabilities needs to be translated to be a node.js module.

Is anyone aware of an existing node.js osr/ogr 'wrapper' (or for generic javascript)?


texts: Scaling DA UpdateCursor to large datasets?

texts:
I am currently having what appears to be a scaling issue with arcpy, da update cursors, and large file geodatabases.

I have a piece of code that iterates through every feature in a feature class and does some calculations and manipulation on the data. It works great for smaller data sets, but it is orders of magnitude slower on large ones. I have a simple counter in it and print statements every 1000 records for the purpose of benchmarking it and here are the results:

48k features, each 1000 takes ~0.12 seconds to do
133k features, each 1000 takes ~0.12 seconds to do
2mil features, each 1000 takes ~0.17 seconds to do
48mil features, each 1000 takes ~23 seconds to do (yes, twenty three seconds, no decimal place)

I only let the 48mil run for a few minutes before I killed it but the rest of them were ran until completion and these times are very consistent, be it the first or the last one with only a few hundredths of a second deviation once in a while. Even if they were scaling linearly I would only expect .05 seconds every 2mil records and that would put 48mil somewhere in the 1.3 second area. The results are 20 times that.

The ones with less features are just subsets of the data for testing purposes so it should not be an issue with differences in the data. All of them were created in the same methods, so I do not think there is an issue with differences in the data other than size. It seems to be coming purely from the number of features.

Sorry I don't have the exact code with me at the moment, it is at work, but this is basically it. The code itself works great though, I am mostly wondering if anyone as run into similar issues with large datasets. I was thinking that it could be a memory leak error, but these are all run outside of Arc in the python window and would expect that a memory leak would take time to ramp up (slow down) over time rather than instantly being much slower.

import arcpy
from datetime import datetime

...

i = 0
fromlast = datetime.now()
with arcpy.da.UpdateCursor(fc, fields) as rows:
    for row in rows:
        if i % 1000 == 0:
            now = datetime.now()
            print i, ': ', now - fromlast
            fromlast = now
        ###do stuff here
        rows.updateRow(row)
        i += 1

I only have ArcGIS (ArcInfo license), gdal/ogr2ogr, and python to work with, but I am not set on using the FGDB or da cursors if there is some better way to do it within my limited selection of tools.


texts: ExecuteError: Failed to execute. Parameters are not valid error when SelectLayerByAttribute run in Python 2.7

texts:
I am running a line of code in a Python 2.7 environment using arcpy.SelectLayerByAttribute_management:

fl_haz_selects = arcpy.SelectLayerByAttribute_management(S_Fld_Haz_Ar, "NEW_SELECTION", "\"DFIRM_ID\" = '{}'".format(fipC))

It runs perfectly fine in Python 3.6. Now I am getting the following errors. What am I doing incorrectly here? Note: S_FLD_Haz_Ar is a shapefile. Do I just need to make it into a layer? I am trying to avoid making a physical layer.

  Traceback (most recent call last):   File "", line
  3, in    File "C:\Program Files
  (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\management.py", line 7713, in
  SelectLayerByAttribute
      raise e ExecuteError: Failed to execute. Parameters are not valid. The value cannot be a feature class ERROR 000840: The value is not a
  Raster Layer. ERROR 000840: The value is not a Mosaic Layer. Failed to
  execute (SelectLayerByAttribute).


texts: Using PostGIS to map multiple GPS tracks

texts:
I have a table which contains GPS track points from multiple vehicles.

What I want to do is to produce GeoJSON polygons to overlay each trip on a map over a range of dates.

My GPS table contains:

id - a unique id
device_id - a unique device id,
vehicle_id - a unique vessel id
timecol - the gps timestamp as a unix time stamp
lat - latitude
lon - longitude
location - geom location
distance_from_depot - km from the depot
distance_from_depot_id - a unique depot_id

When the NMEA strings are are imported to my database, I run a trigger which calculates the location, the nearest depot_id, distance from the nearest depot.

To generate the polygons, I currently have a query which I run which delivers the data to my server, which then has to take that data and manipulate it to GeoJSON features  with properties and geometry; a job which should really be done by the database layer.

My current query is contained within a function:

FUNCTION public.get_gps_tracks(p_vehicle_id integer, p_start_timestamp timestamp with time zone, p_end_timestamp timestamp with time zone)
RETURNS json
  LANGUAGE sql
AS $function$
    select 'GPS' as display_group,
        array_agg(distinct entering_depot) as col4,
        min(timecol) as start_time,
        max(timecol) as end_time,
        avg(lat) as data_a,
        avg(lon) as data_b,
    from gps_data
    where
        vehicle_id=p_vehicle_id
        and timecol > p_start_timestamp
        and timecol < p_end_timestamp
    group by (round(timecol/300) sort by timecol;
$function$

I have had some success using a combination of row_to_json, array_to_json and array_agg, but not enough to get me that last 2% of the the problem.

How could I write this to return valid GeoJSON features, one trip per feature like the following GeoJSON .

Perhaps the answer is to use a window function to detect a change in {col4} and group by that?


texts:
I need help on buying 1 x ArcGIS 10 for Server licenses for a production environment. We have three servers with below configurations:

Server 1 - 8 core, 24GB, 64-bit OS, Windows 2008 Server R2 Enterprise
Server 2 - 8 core, 24GB, 64-bit OS, Windows 2008 Server R2 Enterprise
Server 3 - Database server - Linux OS - Oracle 10g - 64 bit with ArcSDE 10

Our plan is to use Server 1 as a main server, where SOM and SOC will be installed and Server 2 as load balancing server (i.e only SOC will be installed).

Our application consumes 1 map service (with 80 layers) and 2 GP services (50 tools for various GIS operation) and 1 Geometry service. NO EDITING FUNCTIONALITY. Lot of saving in SDE database will be done with other business requirements. Roughly 40 concurrent users at the maximum in any given point of time.

We also use some standard ArcGIS GP Tools like the Aggregate Tool, which warrants an enterprise edition of ArcGIS server.

Can you please explain how an ArcGIS Server License (Enterprise Advanced) needs to be procured?

I am looking for optimized licensing for this setup.


texts: QGIS GUI Python Editor Code Running Exception

texts:
import sys
from PyQt4 import QtGui

class MessageBox(QtGui.QWidget):
    def __init__(self, parent=None):
        QtGui.QWidget.__init__(self, parent)
        self.setGeometry(300, 300, 250, 150)
        self.setWindowTitle('message box')
    def closeEvent(self, event):
        reply = QtGui.QMessageBox.question(self, 'Message',"Are you sure to quit?", QtGui.QMessageBox.Yes,QtGui.QMessageBox.No)
        if reply == QtGui.QMessageBox.Yes:
            event.accept()
        else:
            event.ignore()

app = QtGui.QApplication([])
qb = MessageBox()
qb.show()
app.exec_()

Above code executed in python console editor. it will run and open the form on top of qgis it functioning also good but after close the gui from qgis application not responde it will be idle not functioning any operation. Again close the qgis application and reopen it to do the further operation in code.


texts: How to use a list of numbers from a spreadsheet.csv field in the query of arcpy.Select_analysis

texts:
Have an excel spreadsheet in CSV format with a field of APN numbers (8 digit US assessor parcel numbers)

Need to query these in a large layer of parcel polygons using python script and the arcpy.Select_analysis function.

The problem may be the way I created the list object in python.  When I print the list, commas are returned inside parenthesis' with the APN's .

The error message is invalid expression, here is the script I have:

import arcpy
import sys
import os
from arcpy import env
env.workspace = "D:/Project"
env.overwriteOutput = True

APNlist = []
Plist = "D:/Project/Preserve_Parcels.csv" #change format??
cursor = arcpy.da.SearchCursor(Plist, ["APN8"])

for row in cursor:
    APNlist.append(row) #should I change how the list is appended?
print APNlist

infc = "D:/Project/Final.gdb/Parcels"
fieldname = "APN_8"
outfc = "D:/Project/Final.gdb/parcels_sel"
qry = """APN_8 IN {0}""".format(str(tuple(APNlist))) # Wrong syntax?

arcpy.Select_analysis (infc, outfc, qry)

I'm a python/arcpy beginner


texts: Correctly import .tif files in GRASS GIS

texts:
I am trying to import a GeoTiff (.tif) into GRASS GIS. I then want it to show up in the map display.

This is the data I want to use: http://ghsl.jrc.ec.europa.eu/ghs_pop.php

As far as I can tell, this is raster data.
I know the location and the dataset must share the same projection. When I set up the location, I selected "Read projection and datum terms from georeferenced file". GRASS GIS correctly read projection and datum.

Now, I selected the location I just created, and the predefined mapset "PERMANENT". Then I hit "Start GRASS Session".

But now I'm just seeing a purple rectangle. See this screenshot:

I'm not quite sure what's going wrong.

Background: I am to visualize the population growth in south america.
Our college course didn't do anything to prepare me for this so I have to ask such basic questions.


texts: Python script for calculating NDVI in QGIS

texts:
I am looking for developing python codes mainly for working with raster data, in particular satellite imagery from Landsat. I am familiar with GIS and RS and I usually work with QGIS, but I am new in python programming. I tried to write a script to calculate NDVI, such a first training sample, and I found a python code to make NDVI and to run in QGIS but in my case it doesn't work. Can anyone help me to understand why? I am working with the last version of QGIS Essen 2.14 and the code I used is the following:

#Calculate NDVI for QGIS
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

def NDVI (band4,band5,output):

    entries=[ ]
    #define raster 1 (band4)
    raster1=QgsRasterCalculatorEntry()
    raster1.ref='band4@1'
    raster1.raster=band4
    raster1.bandNumber=1
    entries.append(raster1)
    #define raster 2 (band5)
    raster2=QgsRasterCalculatorEntry()
    raster2.ref='band5@1'
    raster2.raster=band5
    raster2.bandNumber=1
    entries.append(raster2)
    #NDVI Processing
    calc=QgsRasterCalculator('(band5@1+band4@1)/(band5@1-band4@1)',output,'GTiff',band4.extent(),band4.width(),band4.height(),entries)
    calc.processCalculation()
canvas=qgis.utils.iface.mapCanvas()

layers=canvas.layers()
layers[0].name()

layers[1].name()

NDVI(layers[1],layers[0],'C:\Users\.....')

I would like to use it both to run like a script in the Python console of the QGIS, as well as to made a tool for the for the scripts in the processing toolbox (“create new script”).

Is it possible to save it as temporary file in the script, instead of giving a path?

Indeed the script as I wrote runs, but it doesn't produce any output. But, by putting quotes as suggested by Matte, the script works well and produces results although something wrong is happening because the resulting NDVI is not correct. I need to review it.


texts: Converting a column with coordinates in Pandas to GeoPandas


texts: Rasters in R - merge taking mean, excluding NAs. Overlay and Mosaic both giving funny results

texts:
I have a series of raster images, all derived from LandSat data, which give pixels a value of 1 if they are predicted to contain water, 0 otherwise. Many of these images are from the LandSat 7 ETM+ SLC-off, which means they contain bands of missing data (NAs). Note the exact position of those bands varies from image to image. What I want to do is, for each pixel, take the mean value, simply excluding NAs.

I have tried the following (where water_rasters is a list containing the individual images described above):

water_mean <- do.call(overlay, c(water_rasters, fun = mean, na.rm = T)

and

water_mean <- do.call(mosaic, c(water_rasters, fun = mean, na.rm = T)

When using the overlay function, it appears that all pixels where any of my input layers were NA are NA in the output, so the water_mean raster is almost entirely NAs.

When using the mosaic function, I'm not quite sure what it does but the final image is certainly not the mean - almost every pixel is either 0 or 1, with no intermediate values (where I would expect quite a few intermediate values, and I have even checked individual pixels from the inputs and some should definitely be intermediate). However, it does seem to cover the whole scene, and does not have bands of NA data. As I said, I can't tell exactly what mosaic is doing - the water_mean raster certainly is generated from the input data (I can make out rivers, lakes etc) but it's not the mean of all of them or even the mean of a subset of them.

Thoughts? Workarounds? (Would be happy to use QGIS instead, but want the solution to be easily scaled up and repeatable.)

EDIT
As per request, simple example:

raster1 <- raster(matrix(c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1),5,5))
raster2 <- raster(matrix(c(NA,NA,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1),5,5))
raster3 <- raster(matrix(c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0),5,5))
raster4 <- raster(matrix(c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1),5,5))
raster5 <- raster(matrix(c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1),5,5))

extent(raster1) <- c(36,37,-3,-2)
extent(raster2) <- c(36,37,-3,-2)
extent(raster3) <- c(36,37,-3,-2)
extent(raster4) <- c(36,37,-3,-2)
extent(raster5) <- c(36,37,-3,-2)

projection(raster1) <- CRS("+proj=longlat +datum=WGS84")
projection(raster2) <- CRS("+proj=longlat +datum=WGS84")
projection(raster3) <- CRS("+proj=longlat +datum=WGS84")
projection(raster4) <- CRS("+proj=longlat +datum=WGS84")
projection(raster5) <- CRS("+proj=longlat +datum=WGS84")

water_rasters<-list(raster1, raster2, raster3, raster4,raster5)

water_mean_overlay <- do.call(overlay, c(water_rasters, fun = mean, na.rm = T))

water_mean_mosaic <- do.call(mosaic, c(water_rasters, fun = mean, na.rm = T))

In this example, water_mean_mosaic IS the desired output raster, where as overlay doesn't work (has NAs in the two top-left cells, where raster2 is NA).

So I guess the problem isn't my code per se... Could the problem with mosaic be due to large numbers of large files? I'm talking about ~100 LandSat images, each 7011x7611 pixels...


texts: Add layer keywords to WMTS GetCapabilities response using GeoServer?

texts:
In the past I have been entering keywords in "Edit Layer" section of GeoServer. These keywords have shown up on a per layer basis in the "KeywordList" section of a WMS GetCapabilities response, which has allowed me to create a layer tree for my map application. I have recently moved over to using WMTS GetCapabilities requests, however these responses do not show the "KeywordList" section. Is there a way to get this metadata to show up in WMTS GetCapabilities responses?

I have seen there is a "INSPIRE extension" for GeoServer which says it can expand WMTS metadata requests, however after installing it successfully, I have been unable to get the "KeywordList" section to show up. Does anyone have any suggestions?


texts: Using QGIS 2.14.2 after installing qgis2web, the app doesn't appear in the plugin tab, even though all the files were successfully downloaded and verified in the python plugin directory. Other plugins I downloaded have installed just fine. Tried uninstall and install again, tried reinstalling, tried quitting and restarting qgis, still no app.


texts: Privileges, Register as Version commands are greyed out when right-clicking a feature class in an ArcSDE geodatabase

texts:
When right-clicking a feature class in an ArcSDE geodatabase from ArcGIS Desktop, the 'Privileges', ‘Register As Version’ commands are grayed out even though the user is connected as the data owner, and (Advanced) ArcGIS Desktop license is being used.

In my case the data owner is the ‘SDE’ user, I restored this database from the backup to a new SQL server instance on a test machine, The SDE user is able to connect to the database and view the feature classes, now I am trying to assign privileges to users and register the feature classes as versioned and I found out that the Privileges, Register as Version commands are greyed!

I made sure that I am using the SDE user who owns the data and an Advanced license of ArcMAP 10.5.0 https://support.esri.com/en/technical-article/000009946

I tried to use change privileges (data management) tool but I got this error!

ERROR 999999: Error executing function.
DBMS table not found [42000:[Microsoft][ODBC Driver 13 for SQL Server][SQL Server]Cannot find the object 'CanLocations_evw', because it does not exist or you do not have permission.]
Failed to execute (ChangePrivileges).

I checked the view in the database in its there!
What might be the issue here!


texts: Working with WKT data in XML File in SQL Server

texts: Anyone have experience with importing XML files into SQL server 2008 or 2012 that have a WKT element that includes the geometry for each feature?  I can get this to easily work using FME software but wondering if this can't be done directly using SQL Server.  In FME I used the GeometryReplacer transformer. These XML files have all of the WKT data in a single XML element for each feature.


texts: How to properly reorder the raster data from bottom up to top down?

texts:
When creating a raster I set the lower-left x/y values, and asc file has the following "header":

ncols        1731
nrows        1863
xllcorner    17548488.944545675069
yllcorner    -8069349.376163828187
dx           500.162987955558
dy           500.081943659035

Note that the raster data is going from bottom to top.

However, now it seems that I need to reorder the raster data so that my raster data will appear from top to bottom, instead of bottom to top. In the above case, I want the dy to assume negative value.

Is there any built-in function for this?

I could have reorder it manually:

   using (var ds = Gdal.Open(demInputTiff, Access.GA_Update))
    {
        var nCol = ds.RasterXSize;
        var nRow = ds.RasterYSize;
        double[] geoTransform = new double[6];
        ds.GetGeoTransform(geoTransform);

        if (geoTransform[5] < 0)
            return; //no need transformation

        geoTransform[3] += geoTransform[5] * (nRow -1);
        geoTransform[5] *= -1;

        ds.SetGeoTransform(geoTransform);

        Band band = ds.GetRasterBand(1);
        var allValues = new int[nCol * nRow];
        var noData = band.NoDataValue();
        band.ReadRaster(0, 0, nCol, nRow, allValues, nCol, nRow, 0, 0);

        //manual reordering, crude way of iterating through the array
        var allValues2 = new int[nCol * nRow];
        for (int i = 0; i < nRow; i++)
        {
            for (int j = 0; j < nCol; j++)
            {
                allValues2[(nRow - i - 1) * nCol + j] = allValues[i * nCol + j];
            }
        }

        band.WriteRaster(0, 0, nCol, nRow, allValues2, nCol, nRow, 0, 0);

        ds.FlushCache();
    }

But I'm pretty sure that this is not the most efficient way of doing things, both memory and speed wise. There is also an issue of whether .Net array can hold that big the array size or not.

Is there a built-in function for this?


texts: Getting a JSON response with requests and ArcGIS REST API

texts:
I'm performing a batch download from my AGOL account, and using the ESRI documentation export item I can export content to shapefiles.  I would like to then download the newly created shapefiles.  In order to accomplish this I have to be able to have the ID number of the newly created shapefile.  According to the documentation it is in the JSON response.

sudo code:

import requests
import arcrest
import json

exporturl = self.root_url +'/users/' + self.username + '/export'
export = requests.post(exporturl, params=FsParams)
export_response = export.json()

What I have tried:

response = export_response['exportItemId']
print response

Error:  KeyError: 'exportItemId'

printing response = export_response returns the params from the post

Using the requests library can I access the JSON response?  If this is not possible with requests what method will work?


texts: Does spatial indexing speed batch Point in Polygon

texts:
I have two tables, one with 66 thousand points and other with 5.5 thousand polygons. I'm trying to correlate which point is in which polygon and it is taking a long time. I want to know if spatial indexing should improve the query speed.

The polygon layer is already simplified and I've already created the spatial index for both the coordinates layer and the polygons layer.

I've tried this on PostGIS with SQL queries and in QGIS with it's Point in Polygon algorithm. Both of them take several minutes to process.

Does spatial indexing really help with this kind of problem, or just the nearest neighbors or intersections kinds of problems?


texts: Modelize a road networks out of GPS data

texts:
I have a lot of GPS data from vehicles driving around on a private property which is not mapped anywhere.

I would like to model the road network using this GPS data ; That would allow me to check where those vehicles are live by snapping their position back on the road we've just modeled.

I have a big time finding resources online on how this can be done in python for instance. I would think that this is a problem that has been "solved" already and that libraries would be available. Can anyone drive me through the available resources?


texts: Calculating water body area in Google earth engine using Landsat images?

texts:
I am attempting to calculate the area of water bodies  in my classified map using landsat images in Google earth engine. from youtube I found a tutorial on GEE which showed me a method of calculating an area. However I struggled to apply it because I don't understand why my code is showing error.

Can someone guide me through how to calculate area, whether using the same method, or another?

Here's the link to my map:

https://code.earthengine.google.com/ff99d299831a9768fc70cd8fec15807b

the code is:

var l5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR"),
    table =
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[90.95591178825589, 23.73559234169589],
          [90.95591178825589, 23.528001478587907],
          [91.1265431969473, 23.528001478587907],
          [91.1265431969473, 23.73559234169589]]], null, false);

var cloudMaskL457 = function(image) {

  var qa = image.select('pixel_qa');

// If the cloud bit (5) is set and the cloud confidence (7) is high

// or the cloud shadow bit is set (3), then it's a bad pixel.

  var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));

// Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2);
};

l5 = l5.filterDate('1989-09-01', '1990-03-30')
      .filterBounds(table)
      .filterMetadata('CLOUD_COVER_LAND', 'less_than', 10)
      .map(cloudMaskL457);

Map.centerObject(table, 12);
print('l1: ', l5);

var visParams = {bands: ['B3','B2', 'B1'], min: 600, max: 2500};

//['B3_median','B2_median', 'B1_median'], min: 600, max: 2500};

var visParamsFalse = {bands: ['B4','B3', 'B2'], min: 600, max: 2500};

var median = l5
//.reduce(ee.Reducer.median())
      .median()
      .clip(table);

print(median)

//Map.addLayer(median, visParams, 'median');
Map.addLayer(median, visParamsFalse, '1990False');

//Map.addLayer(awei1, {min: -1, max: 1, palette: ['FF0000', '0000FF']},'1990l1');

var awei2 = median.expression(
    'blue + 2.5 * green - 1.5 * (NIR + MIR) - 0.25 * SWIR', {
      'blue': median.select('B1'),
      'green': median.select('B2'),
      'NIR': median.select('B4'),
      'MIR': median.select('B5'),
      'SWIR': median.select('B7')
}).rename('awei2');

print(awei2)

Map.addLayer(awei2, {min: 1, max: 2, palette: ['FF0000', '0000FF']},'classified');

var area_pxa = awei2.multiply(ee.Image.pixelArea())
                    .reduceRegion(ee.Reducer.sum(),table,30,null,null,false,1e13)
                    .get('constant')

 area_pxa = ee.Number(area_pxa).divide(1e6)
 print ('Area using ee.Image.pixelArea (km²)', area_pxa)


texts: Interpolating values within certain distance of points using ArcGIS Spatial Analyst?

texts:
Let's say I have a set of points representing bridges.  I want to create a kind of "heat" map that would highlight areas within my study area that has many bridges in good or bad condition.  I have a rating attribute for each bridge point.

I've tried the spline interpolation method (in ArcGIS spatial analyst) which seems to work okay but it interpolates values for the entire area.  I only want raster values for areas that within a certain distance of a bridge point (say one mile for example).  How would I go about doing this?


texts:
When setting up a queryTask, what's the difference between queryTask.execute() and queryTask.executeForIds()?

In ESRI's JS API Reference, they have the exact same descriptions:

  execute(parameters,callback?,errback?)

  Executes a Query against an ArcGIS Server map layer. The result is
  returned as a FeatureSet. If the query is successful, the
  user-specified callback function is invoked with the result. A
  FeatureSet contains an array of Graphic features, which can be added
  to the map using Map.graphics.add(). This array will not be populated
  if no results are found.

  executeForIds(parameters,callback?,errback?)

  Executes a Query against an ArcGIS Server map layer. The result is
  returned as a FeatureSet. If the query is successful, the
  user-specified callback function is invoked with the result. A
  FeatureSet contains an array of Graphic features, which can be added
  to the map using Map.graphics.add(). This array will not be populated
  if no results are found.

The examples are slightly different, but I can't tell how these methods differ.


texts: GRID-file: get max extent for all z-values above given value

texts:
I have a GRID-file with a quite large extent.
For an automatic workflow I need the maximum extent/bounding-box for all z-values that are over a certain range.
For example: Get the maximum extent for all z-values that are greater or equal than 0.05

Is there a possibility to do that with Python, GDAL or GMT-tools (generic mapping tools), so that it can run automatically?


texts:
I've been attempting to use the qgis module (QGIS version 3.14) in standalone scripting on a Windows 10 machine. However, I have found myself unable to import the qgis module. Code and message below:
import qgis

returns
No module named 'qgis'

How do I either install qgis or find the appropriate python installation on my computer?
I've seen other versions of this question, however, they all seem to be for MacOS or Linux or I've attempted the provided solutions and had no results.
My attempts have generally been within the python window itself, not cmd or qgis console. And I've tried every installed version of python on my machine. even the one within the QGIS install folder.
EDIT: I've discovered within the bin of my QGIS install, python-qgis.bat. running this opens a python window that allows the import qgis to work.
BatchFile code
@echo off
call "%~dp0\o4w_env.bat"
call qt5_env.bat
call py3_env.bat
@echo off
path %OSGEO4W_ROOT%\apps\qgis\bin;%PATH%
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT:\=/%/apps/qgis
set GDAL_FILENAME_IS_UTF8=YES
rem Set VSI cache to be used as buffer, see #6448
set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python;%PYTHONPATH%
"%PYTHONHOME%\python" %*

My question now is, I typically run my python scripts via a batch file. Will I need to include this at the start of my own batch files. I don't entirely understand what this file is doing.


texts:
I have a peculiar rendering of a small table in PostGis of polygones (imported from a shape file) from Geoserver (layer preview with OpenLayers)

The object are just cloned ... Would you give me a clue on what is going on ?

Here is the output of the same table using QGis


texts: Connecting vector label size with map-units?

texts:
I have a working application of labels to a point array - but the size of the text is not relative to the map scale.  This decreases the readability of labels when panning and zooming.

To resolve~

How would one implement a text-size based on MapUnits?

    field = str(csvGridFieldNames[index])
    csvGrid = ftools_utils.getMapLayerByName(unicode('mfLayer1_Grid'))#str(activeLayer)))
    textSize = 7
    valLabel = QgsPalLayerSettings()
    valLabel.readFromLayer(csvGrid)
    valLabel.enabled = True
    valLabel.fieldName = field
    valLabel.placement= QgsPalLayerSettings.OverPoint
    valLabel.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'%f' %(textSize),'')
    valLabel.writeToLayer(csvGrid)
    iface.legendInterface().refreshLayerSymbology(csvGrid)
    iface.mapCanvas().refresh()

This question appears related, but having some trouble understanding the doc's...

"...lables 'Data Defined' section in pyqgis"

...and FWIW:

A screen shot of the Qgis tool, with the two (needed) PyQgis settable parameters highlighted in green:

Updates following underdark's advice:

Input syntax errors (on my part) are not raising error flags - but are also not implementing MapUnit based font/label size...

    valLabel.SizeUnit = QgsPalLayerSettings.MapUnits
    valLabel.setDataDefinedProperty(QgsPalLayerSettings.FontSizeUnit,True,True,':%f' %(textSize),'')

LAST UPDATE [SOLVED]
(final working code below):

    textSize = (abs(scale[0]-scale[1]))/4
    valLabel = QgsPalLayerSettings()
    valLabel.readFromLayer(csvGrid)
    valLabel.enabled = True
    valLabel.fontSizeInMapUnits = True
    valLabel.fieldName = field
    valLabel.placement = QgsPalLayerSettings.OverPoint
    valLabel.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'%f' %(textSize),'')
    valLabel.writeToLayer(csvGrid)
    iface.legendInterface().refreshLayerSymbology(csvGrid)
    iface.mapCanvas().refresh()


texts: How to design a diagram/framework to exhibit tool's underlaying structure

texts:
Possible Duplicate:
  How to generate schemas of GIS hardware / software / workflows?

Can anyone give me some examples of framework diagrams that showcase/exhibit a tool, web-mapping application or similar that use a smorgasbord of interconnected elements, libraries, functions and tools. This is especially useful for communicating how your tool works when speaking publicly (per se) about the tool itself.

I am in the process of creating a poster for an academic conference (GeoIDE Geospatial Conference) and need some help putting my idea onto paper. Note that I may be using the incorrect terminology... feel free to edit!

This is a rudamentary proposal for what I want to display:

             DECISION   ALTERNATIVES
                    ↑   ↑
                     MAP <------------- ↰
                     ↑↑↑                 |
MAPSERVER ----◇-> OPENLAYERS      DECISION MAKER
     JQUERY --↥------------------------ ↵

Is there a best practise to this type of diagram/visual aid? What have you come across that helped you design something similar? Any help is appreciated. I'm about to break out the crayons and construction paper!


texts:
I am trying to plot a bunch of areas on a Google map. The locational elements of the data seem to be in St_area and St_length fields. I can't find much about this format. Is anyone familiar with how I can convert this to a set of coordinates that I can then put on a Google map?

Data here:

https://geoportal.statistics.gov.uk/datasets/da831f80764346889837c72508f046fa_3


texts: using MCA in GIS for suitability analysis

texts:
I have an essay about using MCD methods to solve suitability analysis.
I must write a scenario about fuzzy membership, WLC and AHP.

I feel that I’m confused about when and where to use each method. I read many papers that used the 3 methods in the same problem.
For example in landfill site selection,

first the fuzzy function used for data standardization (0,1).

Here I have a question: after determining the criteria (for example the distance from roads ) how to calculate the fuzzy function)?

Then WLC was used to give the different variables weights>

I need to know how weights are calculated?

Then AHP used to choose to best alternative.

In this case where Boolean can be used and how?
Is it possible to use one method only to solve such problem and get a good result?
Also I would like to know advantages and disadvantages of each method.


texts: Map with style (color, border) depending on property value

texts:
Using Leaflet combined with D3, I'm trying to apply a style (e.g. colorize or draw a border) to units that share a same property value.

Further explanations

Context: townships in France are grouped in local 'federations' ("communauté de communes" used as frameworks to carry local tasks: collecting garbage, managing schools...). It is usually admitted that these federations are formed according to a concept called "bassin de vie" (https://fr.wikipedia.org/wiki/Bassin_de_vie) which is the smallest territorial girth where people live. However, these federations may also be formed according to the global income of their townships's inhabitants, and this is what I'm trying to visualise.

The data I'm using is shown below, and contains the ID of the federation every township is part of.

How can I draw the external borders of the townships being part of the same federation?

The data structure I use:

{"features": [{
    "type": "Feature",
    "geometry": {
        "type": "Polygon",
        "coordinates": [
            [
                [x1, y1],
                [x2, y2], ...
                [x3, y3]]]
    },
    "properties": {
        "zip_code": "xxx",
        "township_name": "xxx",
        "township_income": xxx,
        "township_population": xxx,
        "federation_code": xxx
    }
},


texts: IProximityOperator.ReturnDistance returning incorrect result?

texts:
I am trying to implement a search of nearest features in ArcObjects.NET (don't suggest existing tools like Near) and I use IProximityOperator.ReturnDistance to get distances between features and to select the nearest one.

However, I've found that sometimes this method returns incorrect results.

Here is test map with one point layer (in general type of geometry is arbitrary):

Test Data.

Look at the feature with OID = 18. You can see that the nearest point is a feature with OID = 40 but the result of IProximityOperator.ReturnDistance for this feature is greater than for feature with OID = 22.

What am I doing wrong?

Test code:

var f18 = featureClass.GetFeature(18);
var f22 = featureClass.GetFeature(22);
var f40 = featureClass.GetFeature(40);

var f18AsProximityOp = (IProximityOperator)f18.Shape;

var d18_22 = f18AsProximityOp.ReturnDistance(f22.Shape);
var d18_40 = f18AsProximityOp.ReturnDistance(f40.Shape);

But d18_22 < d18_40. I use ArcGIS 10.2, .NET 3.5, Windows 7.


texts:
I have tried to create  map of south west pacific and have got difficulties in interpolating data using different crs.
the interpolation is broken and have a line in the 180 degrees marking two different sectors of the map
i am trying to fix this as it seems the map is divided into two different sectors
Can some expert assist in this?


texts: How to create a straight line from start and end vertices of a line using Qgis?

texts:
I want to create a straight line from two vertices of a line which are not straight, but I should be getting a straight line joining the start and end vertices of this curve line.

I want to attain the green line


texts: QGIS has duplicated SVG groups in Styling menu

texts:
I'm running QGIS 2.14.8 (latest corporate user version) and SVG groups in the styling menu are duplicated 8 times. I searched for a bug entry but didn't find anything so I thought I'd post here.

Any idea how to clean the duplicates out of the menu?

I checked the instal path and SVG directory only has each folder once.


texts:
I have come code: Is it possible to get coordinates of allready drawn points ? Or get intersection with other feature when drawing?

app.current.interaction = new ol.interaction.Draw({
    features: app.current.featureOverlay.getFeatures(),
    type: featureType.drawType
});
app.current.interaction.on('drawstart', function () {
        app.current.interaction.lastEvent = 'drawstart';
    });


texts: Get daily mean of AOD from MCD19A2 in GEE

texts:
I tried to get daily mean of AOD from MCD19A2 in GEE, but my mean function always returns zero.
// function to add date to image collection from system:time_start
var writeDate = function(image){
  var date = ee.Date(image.get("system:time_start"))
  var day = date.get('day')
  var month = date.get('month')
  var year = date.get('year')
  var image_dated = image.set('day',day)
                          .set('month', month)
                          .set('year', year)
  return image_dated
}
var modis = ee.ImageCollection('MODIS/006/MCD19A2_GRANULES')
                  .select('Optical_Depth_055')
                  .filterDate('2016-05-15','2020-08-31' )
                  .filterBounds(Hanoi)
                  .map(writeDate)
                  ;

var band_viz = {
  min: 0,
  max: 800,
  palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']
};

var first =modis.first()
print ("modis_first",first)
Map.addLayer(first, band_viz, "first")

var date = first.get("system:time_start")
// print("date",date)

var date_easyToRead = (ee.Date(date))
print("date_easyToRead",date_easyToRead)

var day = date_easyToRead.get('day')
print('day',day)

var month = date_easyToRead.get('month')
print('month',month)

var year = date_easyToRead.get('year')
print('year',year)

// var first_dated = first.set('month', month)
//         .set('year', year)
var first_dated = writeDate(first)

print('first_dated',first_dated)

var mean = modis.mean()

print("mean", mean)
var meanWithProperties = mean.copyProperties(modis.first(),["system:time_start"])
print("meanWithProperties", meanWithProperties)

// Map.addLayer(ee.Image(byMonth.first()));
var days = ee.List.sequence(15,18);
var months = ee.List.sequence(5,5);
var years = ee.List.sequence(2016, 2016);

// // var a = months.first()
// print("a",months)

var byMonthYear = ee.ImageCollection.fromImages(
  years.map(function(y) {
    return months.map(function (m) {
      return days.map(function (d) {
      return modis
        .filter(ee.Filter.calendarRange(y, y,  'year'))
        .filter(ee.Filter.calendarRange(m, m,  'month'))
        .filter(ee.Filter.calendarRange(d, d, 'day_of_month'))
        .mean()
        .set('day', d)
        .set('month', m)
        .set('year', y)
        .clip(Hanoi)
        // .set("year", (y+m).getInfo())
        ;
  });
});
}).flatten());
print ("byMonthYear", byMonthYear)
Map.addLayer(byMonthYear.first(),band_viz, 'dailyComp')

var reprojected2 = byMonthYear.toBands();
var areaofinterest = reprojected2.clip(Hanoi);
print("areaofinterest",areaofinterest)

Export.image.toDrive({
  image: areaofinterest,
  description: 'MCD19dailyTest5',
  region: Hanoi,
  scale: 1000
});


texts: Interpolation (TIN) from Python Console in QGIS

texts:
I am referring to the most voted answer by root676 in Calling interpolation plugin from Python console of QGIS.

root676 ends his/her answer with:

  Keep in mind that the QGIS-API is currently rewritten to version 3.0
  and the used interpolation-classes are moved from qgis.analysis to
  qgis.core! This will have a huge impact on the functionality of this
  script so that it must be rewritten for version 3.0!

Now I want to use this code to make a TIN interpolation in QGIS 3. However, I could not find any post with code that is adapted to the QGIS 3 version.

Can anyone help me out?


texts: Importing XML workspace file without losing z and m parameters?

texts:
I am attempting to import an XML Workspace document (database schema) from a file geodatabase into an SDE geodatabase sitting on top of Postgres. The steps I am taking are:

In ArcCatalog - right-click on the file geodatabase --> Export --> XML Workspace document (schema only). Save the XML file.
Inspect the XML file to make sure it makes sense. Below on the left the ZOrigin, ZScale, MOrigin, and MScale values all make sense (either 10000 or -100000).
Right-click on my SDE geodatabase connection as the schema owner --> Import --> XML Workspace document.
I then repeat step 1 on the SDE geodatabase (new geodatabase with no data or schema loaded) so that I can compare the XML schemas.

The XML snippet on the right shows my SDE schema. You will notice that the ZOrigin, ZScale, MOrigin, and MScale values are all now either 0 or 1. For some reason the import process is not preserving these values.

I have tried to simply copy and paste the schemas in ArcCatalog and ArcMap, an using the Distributed Geodatabase import tool, but with no success. When I copy or import the XML schema into a new file geodatabase it works just fine. I have also tried creating a new geodatabase user, but this doesn't produce any different results.

This should be a straight-forward export/import process. Am I missing a step or setting somewhere?


texts: GeoExt (ExtJS + OpenLayers) not displaying Google Map

texts:
I've just started exploring GeoExt.  It displays all the map controls, Google Logo and "Terms of Use". But the map is not rendered. I'm using firebug, and no JavaScript errors are generated. I'm using ext 3.4, OpenLayers 2.11, and GeoExt 1.1.

    <link rel="stylesheet" type="text/css" href="libs/ext/resources/css/ext-all.css">
    <script type="text/javascript" src="libs/ext/adapter/ext/ext-base.js"></script>
    <script type="text/javascript" src="libs/ext/ext-all.js"></script>
    <script src="libs/openlayers/lib/OpenLayers.js"></script>
    <script type="text/javascript" src="libs/geoext/lib/GeoExt.js"></script>
    <script src="http://maps.google.com/maps/api/js?sensor=false"
            type="text/javascript"></script>

    <script type="text/javascript">

    Ext.BLANK_IMAGE_URL = "ext/resources/images/default/s.gif";
    var app, items = [], controls = [];

    Ext.onReady(function() {
        app = new Ext.Viewport({
            layout: "border",
            items: items
        });
    });

    items.push({
        xtype: "gx_mappanel",
        ref: "mapPanel",
        region: "center",
        map: {
            numZoomLevels: 20,
            controls: controls
        },
        extent: OpenLayers.Bounds.fromArray([
            -122.911, 42.291,
            -122.787,42.398
        ]),
        layers: [new OpenLayers.Layer.Google(
                'Google Streets',
                {isBaseLayer: true}
        )]
    });
    controls.push(
        new OpenLayers.Control.Navigation(),
        new OpenLayers.Control.Attribution(),
        new OpenLayers.Control.PanPanel(),
        new OpenLayers.Control.ZoomPanel()
    );

    </script>
</head>
<body>
</body>

Any assistance will be appreciated.


texts: Creating line from points using PostGIS?

texts:
I have a PostGIS table with position data of numerous vehicles and want to create lines from these points.

Data can be selected by vehicle id and ordered by timestamp but how to create lines from the result?

What I basically need is one line segment from point 1 to point 2, finalize the line and then again from point 2 to point 3. Of course all that under consideration of the vehicle id.

The latter one is needed because I want to calculate the cruise direction and speed of the vehicle from one point to the next.


texts: Update graphic on Graphic layer with new Geo point (new current location)

texts:
We have migrated ArcGIS library from 10x to 100.9.0.
In our App we have the functionality of change mock location using long press. Below is the code snippet of 10x which is working fine.
Graphic graphic = mGraphicsLayer.getGraphic(marker);
if (graphic != null) {
Point point = LocationUtils.convertToMapPoint(position, getSpatialReference());
final Point geometry = (Point) graphic.getGeometry();
geometry.setXY(point.getX(), point.getY());

mGraphicsLayer.updateGraphic(marker, geometry);
mGraphicsLayer.bringToFront(marker);
}

After Migration we does not found any method like update graphic and change location of marker which is already drawn on Graphic layer. We have more than 1000 marker on Map so also not feasible to clear and redraw all again.
I tried so much but there are no such documentation help. So we are stuck now and do not know how to migrate above code with new library.


texts: 'Snapping to segments' does not work in QGIS 3.10 and 3.4

texts:
QGIS 3.x does not let me snap to vertice and segment anymore, only vertices.

I changed the snapping options several times in the toolbar and in the options menu. Made sure the correct layers are activated. Restartet everything.

I even uninstalled QGIS 3.4 deleted all app data and installed 3.10 instead but the problem prevails.

I can not bring QGIS to snap to segment or follow an existing shape.

Does anyone have a solution?

I am familiar with the settings and choosing the right ones is unfortunately not the problem. Also zooming closer, or starting a new layer did not help. Neither did importing the layers into a new project.

Using standard CRS WGS84

thank you for your answer but i am familiar with the settings and choosing the right ones is unfortunately not the problem.
Also zooming closer, or starting a new layer did not help.
Neither did importing the layers into a new project.

Using standard CRS WGS84


texts: Export images in google earth engine with date name

texts:
I would like to export the NDVI with the date in the filename. Now I can only export them with a list of numbers, e.g. example NDVI_1, NDVI_2 (as in the code that I show below) ... and I want it so NDVI_18042018.

// Load a FeatureCollection
var violada = ee.FeatureCollection("ft:1KmH70D7VKdiWocelf0RtbL_kQhGk0LGkQQ7O-ceG");

// Load a FeatureCollection Sentinel 2.
var s2 = ee.ImageCollection('COPERNICUS/S2')
      .filterBounds(violada)
      .filterDate('2016-10-01', '2017-09-30');
print(s2);

// Function to mask clouds using the Sentinel-2 QA band..
function maskS2clouds(image) {
         var qa = image.select('QA60');
 // Bits 10 and 11 are clouds and cirrus, respectively.
         var cloudBitMask = ee.Number(2).pow(10).int();
         var cirrusBitMask = ee.Number(2).pow(11).int();
 // Both flags should be set to zero, indicating clear conditions.
         var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
         qa.bitwiseAnd(cirrusBitMask).eq(0));
// Return the masked and scaled data.
   return image.updateMask(mask).divide(10000);
}
   var imagen=s2.map(maskS2clouds);

// visualize the first image in the collection, pre- and post- mask
var visParams = {bands: ['B8','B4','B3']}
Map.addLayer(ee.Image(imagen.first()), visParams, 'clouds masked', false)
Map.addLayer(ee.Image(s2.first()), visParams, 'original', false)

//add ndwi
var addNDVI = function(image) {
  return image.addBands(image.normalizedDifference(['B8','B4']).rename('NDVI'));
};
var imagen2=imagen.map(addNDVI);
var NDVI = imagen2.select(['NDVI']);
print(ee.Image(NDVI.first()));
var composite = NDVI.qualityMosaic('NDVI').clip(violada);
print(composite);

// Visualize NDVI
var ndviPalette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163',   '99B718',
  '74A901', '66A000', '529400', '3E8601', '207401', '056201',
  '004C00', '023B01', '012E01', '011D01', '011301'];
Map.addLayer(composite.select('NDVI'),
             {min:-1, max: 1, palette: ndviPalette}, 'ndvi');
//List NDVI
var list=NDVI.toList(150);
for (var i=0;i<150;i++){
     var image=ee.Image(list.get(i));
     var name=ee.String('NDVI_')
                .cat(ee.String(ee.Number(i)))
                .getInfo();
     print(name);
// Export NDVI
     Export.image.toDrive({
           image: image,
           description:name,
           scale: 10,
           region:violada,
           crs : 'EPSG:32630'
           });
}
Map.setCenter(-0.6, 42, 12);

Google earth engine code is: https://code.earthengine.google.com/7101729aaaf2f8be4a0114561b9a673e


texts:
I have multiple large rasters. Their 5 percentile and 95 percentile are needed to be calculated. After loading the rasters to the ArcMap, I am using the following codes -

import numpy
import arcpy

rasters = ["raster1.tif", "raster2.tif", "raster3.tif", "raster4.tif"]

for i in rasters:
    array = arcpy.RasterToNumPyArray(i,nodata_to_value=10000000000)
    marray = numpy.ma.masked_values(array,10000000000)
    del array
    percentile5 = numpy.percentile(marray.compressed(),(5))
    percentile95 = numpy.percentile(marray.compressed(),(95))
    del marray
    print percentile5 , percentile95

For the first two rasters, it provides the percentile values. For the third raster, it is showing memory errors as provided below.

Traceback (most recent call last):
  File "E:\test\Stat\Inlinear data.py", line 89, in <module>
    marray = numpy.ma.masked_values(array,10000000000)
  File "C:\Python27\ArcGIS10.5\lib\site-packages\numpy\ma\core.py", line 2200, in masked_values
    condition = umath.less_equal(mabs(xnew - value), atol + rtol * mabs(value))
MemoryError

Failed to execute (IntelinAir).

I am using ArcMap 10.5 in desktop with 32 GB Ram.

Is there any way to get rid of this memory issue?


texts: ArcSDE installation with PostgreSQL (error - 1015)

texts:
I am trying to install ArcSDE for Postgresql on CentOS 5.5 in i686 architecture. I have successfully installed and configured Postgresql. I did SDE pre-installation and as ArcSDE user, run the ArcSDE 10 installation script that went well without any error. During post-installation, while running $SDEHOME/bin/sdesetup, I am getting error 1015. The log file is showing SDE DBMS Error and is unable to install ST_GEOMETRY. I can see that st_geometry.so is created in /sdeexe100/pg841_st_lib/ directory. I am confused if I missed something during installation/configuration process or if I have to copy st_geometry.so file somewhere else for the installation. I have attached my log file.

Any help is highly appreciated!

Mukesh


texts:
Here is a very specific problem :

I have some mercator flight simulator maps I projected with real coordinates. However, due to game simplification, these maps are not oriented correctly and real bearings don't correspond to ingame ones.

In the following map, arrows show the ingame North (0°) and the number indicates the projected real azimuth.

But my flight plan tool need to display ingame bearings to be accurate and useful !

So, is there a way, or a mathematic function allowing me to recalculate ingame bearings from the map projection, when the figure is in the map bounding box ?

I have the WGS84 coordinates of the 4 corners, and maps are always perfectly oriented NSEW ingame.

I can imagine there could be some maths in the solution.


texts: Using Overpass API find ways inside multipolygon relations

texts:
I'm trying to find all buildings in ways AND relations, type=multipolygon.  The code below works but only finds buildings in ways but not in the relation, type=multipolygon for 'landuse=military'.  I know how to output the rel[landuse=military], but not how to feed those relations nodes into the is_in with the ways nodes.

[out:json];
(

way[landuse=military]({{bbox}});

);

node(w);is_in;
area._[landuse=military];
(
     way(area)[building];
     node(w);
);
out;


texts: FME spatial join by PointOnAreaOverlayer transformer rejects Invalid Polygon Geometry Type

texts:
I would like to perform a spatial join of points and polygons in FME Workbench by PointOnAreaOverlayer transformer. The result should be a country polygons with additional attributes from points. I have a workspace:

What I have in rejected features, in attribute table column fme_rejection_code: Invalid Polygon Geometry Type:

Do you know how can I repair or change the Invalid Polygon Geometry Type features? I have tried GeometryValidator transformer - without results - all areas were rejected in PointOnAreaOverlayer transformer.


texts: Recognized, calibrated test points for a rangefinding function?

texts:
I'm trying to find out if there are scientifically accepted calibration points to measure the effectiveness of a geodetic distance finding function. So far I've been using the listed latitude / longitude coordinates of various airports but I'm not sure where exactly in the airports they're measuring (the control towers?) and its hard to judge my error rate with that unknown.

I thought it might be better if there were points that were a known distance apart - such as where the Prime Meridian crosses Greenwich or the compass rose at the Hoover Dam?


texts: Regrid a NetCDF file to match another without multiplication of pixels (CDO or R)

texts:
I have a NetCDF file (file1.nc) displaying rasterized information for a series of point locations. Each pixel with a number is surrounded by no data (see the colored pixels in the image below). I need to remap this file to match the grid of an additional NetCDF (file2.nc - see b/w pixels in the figure below) without changing the information contained (no interpolation of other math required) in file1.nc.

I tried to do that using CDO, using the REMAP function. In order to maintain the information intact, I selected the -remapnn command assigning to the nearest neighbor in file2.nc the value in file1.nc.
cdo -remapnn,file2.nc file1.nc output.nc

Unfortunately, this procedure splits the pixels in file1 when they equally cover pixels in file2 (see below).

How can I prevent this duplication of pixels? is there a method in CDO or R that could prevent this from happening?


texts: Seeking Louisiana PLSS GIS data - Converting (Township Range & Section) to Lat Long coordinates in 2018

texts:
I have seen this question asked before for other US States, but based upon my understanding of Public Land Survey System, each State has its own reference grid. It appears that the US Government once provided a webservice conversion tool called "Geocommunicator", but this is now offline, denigrated, or no longer supported. Is this BLM document of State PLSS geodatabases meant to be a replacement?

To be specific, I need to convert the Township Range Section boundary information for two properties to Lat Long in the State of Louisiana. These properties are in the swamp and are not accessible by road.


texts: Deleting Oracle SDE features classes properly

texts:
I have got a few feature classes in an SDE connection. This is the situation:

a) They are listed in ArcCatalog (10.1), but the open fails. The message says:

'Error opening feature class, DBMS table not found [ORA-04043: object XXX does not exist]'

b) Obviously, the tool Delete_management cannot be used as the geometry does not exist and it fails too.

c) If you go to your Oracle schema the individual tables corresponding to those geometries are not there, they have been deleted I am afraid and the whole thing makes sense. But under the SDE user, if you open the 'LAYERS' table, the table names are there.

I have read this post, but it does not help me:
Delete SDE feature class features with ArcPy

So, in order to clean the database, shall I delete the entry in the LAYERS table and problem solved? what else or what other steps in Oracle are necessary? every time I open ArcCatalog the corrupted layers are listed and I want to remove them completely.

[1]: Delete SDE feature class features with ArcPy


texts: Merging Degree Two nodes in Stream Network using QGIS 3.4

texts:
I have a stream network and I want to convert it in graph like data structure, for which I need to have precise set of FROM and TO vertices.The geometry of stream network is not correct i.e. the stream segments are not split at intersection on stream network.So, I have to merge the degree two nodes in the network so that my network will have only 1 and 3 degree vertices. The label in the figure represents the stream vertex id and degree of the vertex (ID-Degree) . The degree is the total amount of incoming or outgoing edge in the network.I have two shape-files, one is point layer having vertex ID and degree of each vertex and other is stream network having TO and FROM attributes for each edge but with degree 2 vertex as well which is to be eliminated.

I want it to convert as below.

My approach:- For each edge , I will check the degree of the respective vertex ( assuming any starting vertex) and if the degree is not either 1 or 3 (degree greater than 3 is not expected) then the respective edge is merged and the iteration continues else the next vertex pair is considered.

I think there must be a easier way available as a plugin or Postgres/Postgis based query,please point out to any such tool or modified approach to achieve the same.


texts: How to calculate intersection of two arcs?

texts:
Given two arcs (defined by start and end points with latitude and longitude) how to check if there is intersection (cross point) and if there is -- how to calculate it?

Calculate "manually" -- i.e. I am looking for low-level code, like JS, C#, Python or appropriate equations, nothing high-level, like compute_intersection(a,b,c,d).


texts: gdal/python: TIFF not recognised a supported file format

texts:
I'm trying to open a TIFF using GDAL, read it as a array, do some math, and write it back as a TIFF.

However, when I try to gdal.Open(tif) in Python, I get the error:

RuntimeError: /att/gpfsfs/userfs01/ppl/mwooten3/IDS/climate_data_analysis/MODIS_EVI/AverageOutputs/MODIS_EVI_Jan2001-2015_Average_ave.tif' not recognised as a supported file format.

I just did the same thing the same server, with the same versions of gdal/python, but with a different TIFF, and it worked. Also I looked at the TIFF I'm trying to open and it's not corrupt. Also also, triple checked the path and it is in the right spot.

Any idea why this is happening? What am I missing? I'm on a Linux server using GDAL 2.0.1.


texts: use existing function in arcmap use VBA in Arcobjects

texts: I have a UI tool button, I would like to get and use the existing function (select elements tool) in Arcmap when I activate this button, does anybody know that how to call the existing function using Vb arcobjects ??


texts:
I really like the html annotation tool in QGIS. My plan is to show a photo  (web based image// jpg) on each point of my point shapefile.
Unfortunately labels don't work as they do not parse HTML and the tooltips won't last for long (only on hoover, right?) and can't show external images like from a website. Images need to be local.

So is there a way to create a HTMLAnnotation using pyqgis? I mean, I can
see the HTMLAnnotation block in the qgs file:

  <HtmlAnnotationItem hasFeature="0" feature="0" htmlfile="/home/ricckli/Desktop/test.html" vectorLayer="Accessibility20160219224444279">
<AnnotationItem frameHeight="360" canvasPosX="281" mapPosX="1514906.43412938853725791" canvasPosY="195" mapPositionFixed="1" mapPosY="6873695.00431403797119856" frameColorAlpha="255" frameBackgroundColorAlpha="255" frameBorderWidth="1" visible="1" offsetX="50" offsetY="-50" frameWidth="527" frameBackgroundColor="#ffffff" frameColor="#000000">
  <symbol alpha="1" clip_to_extent="1" type="marker" name="marker symbol">
    <layer pass="0" class="SimpleMarker" locked="0">
      <prop k="angle" v="0"/>
      <prop k="color" v="255,0,0,255"/>
      <prop k="horizontal_anchor_point" v="1"/>
      <prop k="name" v="circle"/>
      <prop k="offset" v="0,0"/>
      <prop k="offset_map_unit_scale" v="0,0,0,0,0,0"/>
      <prop k="offset_unit" v="MM"/>
      <prop k="outline_color" v="0,0,0,255"/>
      <prop k="outline_style" v="solid"/>
      <prop k="outline_width" v="0"/>
      <prop k="outline_width_map_unit_scale" v="0,0,0,0,0,0"/>
      <prop k="outline_width_unit" v="MM"/>
      <prop k="scale_method" v="diameter"/>
      <prop k="size" v="2"/>
      <prop k="size_map_unit_scale" v="0,0,0,0,0,0"/>
      <prop k="size_unit" v="MM"/>
      <prop k="vertical_anchor_point" v="1"/>
    </layer>
  </symbol>
</AnnotationItem>

But I don't want to write it for my own for all the points. especially the  canavasPosX and canvasPosy parameter...


texts: Find row with a polygon containing a specific lat/long

texts:
I have a Postgres table that has two columns – address, and geom (which is a multipolygon without any SRID set, but it should be 2263). I'm trying to find all of the rows that contain a given point, which I only have in lat/long form.

The query I want should look, I think, something like this, but I can't tell if I'm writing it wrong, or I simply don't have any polygons that contain the given point. The point is 40.741431, -74.003334 – the center of Google's building in Manhattan.

SELECT t1.address
FROM t1
WHERE ST_Contains(ST_SetSRID(geometry, 2263), ST_GeometryFromText('POINT(-74.003334 40.741431)', 2263));

What am I doing wrong? I suspect it's an SRID thing... (Or am I doing nothing wrong and my data is just bad?)


texts: SAFE FME: How to iterate through Featureclasses in a GDB and write their Filename to an attribute in a separate featureclass

texts:
I have a GDB which contains raster datasets and Vector featureclasses .
Each Featureclass contains only one polygon. The Polygon represents project areas whose name is also the title of the featureclasse e.g "BP_2247"

In the Featureclass itself this name is not written to an attribute.
I therefore want to merge all featureclasses together whilst also writing their filename into a new attribute in the merge file.

My first problem is how to get Safe FME to iterate through a GDB and only use Featureclasses?


texts: Niche MapInfo Question : Digitising options

texts:
I've recently started using mapinfo to digitise on a daily basis, I am usually georefencing old deeds (on top of OS Mastermap) and then digitising them.

The problem is when drawing the polygon the default colour of that line is grey, it would be much better if that line was a bright colour as it can be easy to lose ones place when digitising from black and white deeds.

So my question is, is there an option to change the colour of the lines of polygon before it's finished.

Hope that makes sense.


texts: Clipping a vector shapefile fails to include all data

texts:
I have a vector file with 4 levels of data within it. This started life as mapinfo tab files (7 files, map,mif,mid etc.) These were passed through mapinfo universal translator, reading original crs (epsg27700) and exporting as .shp file (with its colleagues) again forcing epsg27700.

I know this could be achieved with qgis, but I keep finding comments about additional tasks to do to ensure the CRS sticks with the file. Also Noosa defaults to employ OTF correction so I'm never sure if it is really correct or not. This can be styled with a thematic to depict the 4 layers without issue.

However the shape file is to be clipped by another polygon, also .shp, with same CRS. This is a county border outline.

The vector clip tool in Noosa reports 3 problems with missing CRS, but as I have set to ignore issues, It produces a county shaped output but with only one data layer.

I tested other similar files in Noosa, Las Palmas and Coruna that clip easily - so thought not to be clip tool related. Therefore I suspect something faulty with original file, but of course it displays ok before the clipping. Any ideas welcome for gis amateur.

using Qgis Noosa, win10.
Thanks.


texts: DWG Import: How to resolve vector  too long?

texts: I'm attempting to import a schematic layout however in the DWG/DXF import menu, I get the  "vector  too long" error message for ~1400 entities. It's my first time importing a DWG file, so I'm unsure how to proceed.


texts: Using GRASS GIS in Python standalone scripts

texts:
I am trying to find a least cost path script for Python.
I came across r.drain from GRASS GIS.
Is it possible to use GRASS GIS in a Python standalone script?
If so where can I download the bindings?


texts:
NDVI values raster computed by using the band 3 and band 4 reflectance values b4-b3/b4+b3. The output raster is of 32bit type formate although the possible ranges that shows in statistics parameter are from -0.35 to +0.65. Are this is due to abnormal values that have been accounted for no data values. If yes then how can i remove those values from the output raster or make them out during the calculation of NDVI in raster calculator.

@Whuber. commented in a post normalize raster, but cannot understand how to do into my case of ration images. there are 5 images that i took for comparison but how do i normalize them to compare easily


texts: ArcGIS python "add field" and "calculate field" workflow

texts:
I am new in Python and I am creating a Script. I've crafted a functional code, but I'd like to know if there is an easier way to do it, just for curiosity (And learning new ways!).

The purpose of the Script is to create 5 fields in a polygon shape (previously drawn) and to calculate their values.

My code is this:

import arcpy

env.workspace = "C:/Users/...."

#First, create the fieldlist with a loop for. With this I'll create the table
fieldlist = ["Area","Longitude","Latitude","MaxDepth","MinDepth"]
for field in fieldlist:
    arcpy.AddField_management(My_shapefile,field,"dobule")

#And now calculate fields
arcpy.CalculateField_management(My_shapefile,"Area","!shape.area@squarekilometers!","PYTHON_9.3")
arcpy.CalculateField_management(My_shapefile,"Longitude","!shape.extent.XMax!","PYTHON_9.3")
arcpy.CalculateField_management(My_shapefile,"Latitude","!shape.extent.YMax!","PYTHON_9.3")
arcpy.CalculateField_management(My_shapefile,"MaxDepth","!shape.extent.ZMin!","PYTHON_9.3")
arcpy.CalculateField_management(My_shapefile,"MinDepth","!shape.extent.ZMax!","PYTHON_9.3")


texts: Periodically updating OpenStreetMap tile server to keep it up to date?

texts:
I have successfully built my tile server by following this link (https://switch2osm.org/manually-building-a-tile-server-18-04-lts/) and open layers for viewing the tiles.

How do I keep automate the process of updating my tile server periodically from OpenStreetMap server?

I couldn't find any tutorial related to the same.

I read that Mapnik has tools that pulls data from the OpenStreetMap periodically but could not figure out how to get it up and running.
I also tried to do it using the updating section of (https://switch2osm.org/serving-tiles/building-a-tile-server-from-packages/) but it wasn't helpful too.


texts: Which labelling tool/engine in QGIS will be supported in the long term?

texts:
Why does QGIS appear to have two labelling tools?
(Note this is not new and old symbology, but that...) there are two different GUIs, one in the toolbar (which gives the Data Defined Settings option) and then in the Layer Properties a tab which leads to Advanced | Date Defined Placement + others.
I can turn both on and have two indenpendent lables.

I ask, because what will be supported in the long term?

ciao


texts: Creating custom function in QGIS with self intersects

texts:
I'm trying to create some customs functions in QGIS to automate some works I'm doing in my job.
I have some experience using QGIS expression but I'm trying to create some custom functions now.

I have a layer of points with covid cases and I'm agreggating this points suming the cases depending on the scale of map canvas, as you can see in the expression below:

aggregate( @layer_name, 'sum', "susp_atv" , intersects($geometry, buffer(geometry(@parent),0.004 @map_scale ))

Now I'm having trouble trying to transform this in a custom function.
The function below checks if all geometries of a layer intersects with any geometries of another layer.

@qgsfunction(args='auto', group='Custom')
def intersect(layer_name, feature, parent):

    teste = False
    geom = feature.geometry()
    layer = QgsProject.instance.mapLayersByName(layer_name)[0]
    for feat in layer.getFeatures():
        if geom.intersects(feat.geometry()):
            teste=True
            break

    return teste

How can I update this code to perform the self intersect buffer that I need.


texts: Openlayers Animation Examples and Algorithms

texts:
i wonder for a long time if people in this community uses openlayers for animation needs.
we have know this that animating on browser is a little bit exhausting for it but when a map meets with animation, I believe that good things will be out.

Do you know some good websites or algorithms for animation compatible with openlayers?

these are what I know about animation, some codes and examples. pls share your knowledge about animation with openlayers for people who needs this...

The basic javascript code for adding features on map with time delay..

 time = 0;
    for (var i = 0; i < 5; i++) {
        time += 1000;
        setTimeout(function(j) {
            return function() {
                console.log("var is now", j);
            }
        }(i), time);
    }

i2maps

is a geocomputing environment. It provides a flexible programming framework for knowledge discovery from spatio-temporal data and web-oriented visualisation. The Javascript library is used for building the interactive user-interface, and is built on top of OpenLayers.

OpenLayers Timeline

is a simple library to create a time related map with OpenLayers.
It supports cluster strategy (OpenLayers.Strategy.Cluster) and time filtering.
Currently it supports source data in GeoJSON or GeoRSS format.

Timemap.js

is a Javascript library to help use online maps, including Google, OpenLayers, and Bing, with a SIMILE timeline. The library allows you to load one or more datasets in JSON, KML, or GeoRSS onto both a map and a timeline simultaneously. By default, only items in the visible range of the timeline are displayed on the map.

Bike Share Map

Data is updated automatically every 2-20 minutes from the provider's website. Bike usage is simultaneous usage and includes cycle redistribution. Actual total usage across the day may be much higher. Distribution imbalance - the number of cycles that would need to be moved to a different stand, in order for all stands to be the same % full. Higher numbers indicate a more unbalanced distribution, e.g. many bikes in the centre, few on the edge. This website may be subject to interruption.

Data: Provider's website. Map data: OpenStreetMap contributors, including OS Open Data. Hosting: UCL. Powered by Mapnik and OpenLayers.

RaphaëlJS and OpenLayers Integration.

RaphaëlJS is a small JavaScript library that should simplify your work with vector graphics on the web. If you want to create your own specific chart or image crop and rotate widget, for example, you can achieve it simply and easily with this library.

OpenLayers Filter Strategy

This example uses a filter strategy to limit the features that are passed to a layer. Features bound for this layer have a when attribute with date values. A filter strategy is constructed with a between filter that limits the span of dates shown. A simple animation cycles through the domain of the when values, calling setFilter on the strategy with an updated filter.


texts: Calculate volume of a polygon using leaflet

texts:
I have a polygon and a Digital terrain model, I need to find a way to calculate the 3d volume of the area inside the polygon on leaflet.
Is it possible to achieve it using some sort of formula or do I need an external library?


texts: Bulk KML from shapefile convertion with modelbuilder

texts:
I am trying to batch export kml from shapefiles that are all in WGS84. For this I have created a simple model in model builder to iterate through all feature classes in a directory, make a feature layer and convert it to kml. This is how it looks like.

The first kmz comes out OK. However all the next ones get mixed up and this is what it looks like when I import it to GE (the grey lines going round the Earth).

I am using ArcMap 10.3.1


texts: Saving project with layers as one compiled file?

texts:
Is there a way to save a QGIS project as one large file that includes the layers it points to?

I am working with Tecplot quite a bit, which has the option of referencing individual data files or saving all the data as one compiled file. I was wondering if QGIS has a similar capability.

As long as layer locations don't change, I can simply point to the absolute layer location. However, if I want to create a backup disk or hand over the project to a third party, the folder structure may not be preserved.

I found some discussions about QCompile and GeoPackage, which should do what I want, but from what I understand they are still in development stage.


texts: I need help creating a select statement for an attribute table of temperature data. I'm trying to select only those results where the temperature difference between weather stations was 20 degrees or more in the same day. I'm stumped on how to compare the Max_Temp values across 18 different stations. I have an attributes table with fields for Station_Name, Date and Max_Temp and I would like to select only those results where the difference between Max_Temp >= 20. There are 18 Station_Names and 18 Max_Temps for each of 365 Dates.


texts: Dissolving shapefile but retain attribute fields using ArcGIS Desktop?

texts:
Is there any way to dissolve shapefile but keep attributes fields in the same time?

This is how I want to retain fields:

There are two fields I want to conserve in a shapefile. The first column is Net_ID with long type, and the second column is Geology with string type.

I would like to dissolve the shapefile using field NET_ID. Besides, I would like to keep the Geology field, like the below figure.


texts: Check if a map region has been already visited

texts:
I'm building an iOS app and this is my first experience with maps. The issue is pretty simple: the user can move the map with his/her finger. Each time this happens I download some info from a server (info related to the map region just shown) to display markers. What I want is to avoid reloading and reloading again the same info if the user come back to a previously visited region.

Is there a way to check this?

Maybe an algorithm. I think this is a common issue, but I can't seem to find anything on internet.


texts: Deriving output filename from another parameter in Python Script Tool?

texts:
In an ArcGIS Toolbox script tool, how can I get the default value for one parameter to be derived from the value of another parameter?

GPX To Features does this, but I haven't been able to figure out how. There is logic under Tool >> Properties >> Validation to sanitize the filename but nothing about setting a default value. I've copied the validation code verbatim into my tool but that doesn't work.


texts: Create point data from the attributes of other point layers

texts:
I'd like to create a shapefile point layer (for later interpolation) from the attributes of some other point/line layers I have. Basically I've created a natural neighbour grid of stress direction (I'm a geologist), and converted that to a rather dense shapefile point layer, with an attribute describing direction of stress (angle/180) at each point.

I then have a shapefile line layer, which I've used the field calculator to create a 'direction' (angle/360) for each line segment.

I'm interested in the difference between the direction of each line segment and the direction of stress. Ideally I'd like to be able to create a point layer, with a point at each place the stress grid overlaps (of is near to) a line segment, the attribute I want being the difference between line direction and stress direction at each point. The final output I'm after being a natural neighbour interpolation (I can do this bit) of the angle between my lines and the stress direction.


texts:
I'm working on this script in PyQGIS to take points from a layer called LEK_ALL and snap them to the closest feature of a line layer called LEID_filtered which fulfills a certain request.

In order to do so, I created the function snapPointToLine. This is the full code:

from qgis.PyQt.QtGui import *
from qgis.PyQt.QtWidgets import *
from qgis.PyQt.QtCore import *
from qgis.core import *
import os

def snapPointToLine(feature, lineLayer, request, AttrID):
        '''
        QgsFeature (point), QgsLayer (lines), QgsFeatureRequest, list -> QgsFeature
        Takes the feature feature and snaps it to the colest point of the feature in layer lineLayer that fullfills request.
        Returns the new feature with the same attributes as the old ones, plus the attributes of the snapped-to feature with indices in AttrID.
        '''
        pt = feature.geometry().asPoint()
        lineIter = lineLayer.getFeatures(request) #iterator of the features in lineLayer fulfilling the request
        spIndex = QgsSpatialIndex(lineIter) #create spatial index from selected feature
        nearestIds = spIndex.nearestNeighbor(pt,1) #list of 1 closest neighbour
        featureId = nearestIds[0] #item 0 (and only) of the list is the id
        nearestIterator = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId)) #generate iterator containing only the closest feature
        nearFeature = QgsFeature() #generate an empty feature
        nearestIterator.nextFeature(nearFeature) #copy information from nearestIterator to nearFeature
        closeSegResult = nearFeature.geometry().closestSegmentWithContext(pt) #results in a description of the closest segment between the nearFeature and point
        closePoint = closeSegResult[1] #closest point on geometry to LEK
        #generate feature in new location
        newFeature = QgsFeature()
        newFeature.setGeometry(QgsGeometry.fromPoint(QgsPoint(closePoint[0],closePoint[1])))
        newFeature.setAttributes(feature.attributes() + [nearFeature.attributes()[i] for i in AttrID])
        return newFeature

registry = QgsMapLayerRegistry.instance()
LEK = registry.mapLayersByName('LEK_ALL')[0]
LEID = registry.mapLayersByName('LEID_filtered')[0]

LEKIter =  LEK.getFeatures()
Date_id = LEK.fieldNameIndex('Date') #gets id of field "Date"
LEID_id = LEID.fieldNameIndex('AssetNum') #gets id for field "AssetNum"

newFeat_list = [] #generate new empt list to save the new features

for feature in LEKIter:
        LEKDate = feature.attributes()[Date_id] #get failure date
        y,m,d = LEKDate.year(),LEKDate.month(),LEKDate.day() #split it into components
        request_str = u'"InstYear" <= to_date(\'%d-%02d-%02d\') AND "RemDate*" >= to_date(\'%d-%02d-%02d\')' % (y,m,d,y,m,d) #generate query pt.1
        request = QgsFeatureRequest().setFilterExpression(request_str) #generate query pt.2
        newFeature = snapPointToLine(feature, LEID, request, [LEID_id])
        newFeat_list.append(newFeature)

When I run it I get the following error:

NameError: global name 'QgsSpatialIndex' is not defined

Associated to the line

spIndex = QgsSpatialIndex(lineIter) #create spatial index from selected feature

Any suggestions?


texts: Group members add their location to OSM map?

texts: We have a Facebook group with about 100,000 members. How can we create an OpenStreetMap map where members who have the URL to the map can add their location with a marker (only one marker per member), but not edit the map settings or change markers that other members have added?


texts:
How to write OpenLayers code to create clusters from MVT (using 100k points)?

Because MVT file format saves disk space (90 % less than GeoJSON) and requires only 50 % RAM of browser.

Class hierarchy of mbtiles as vectortile:

  ol.Map / layers: ol.layer.VectorTile / source: ol.source.VectorTile /
  format: ol.format.MVT

Class hierarchy of GeoJSON as cluster:

  ol.Map / layers: ol.layer.Vector / source: ol.source.Cluster / source:
  ol.source.Vector / format: ol.format.GeoJSON


texts: Why are text values ending up in wingdings or something similarly useless?

texts: I have dissolved a layer in QGIS 2.8. The resulting layer has all the column labels and text fields changed to wingdings or something similarly useless. When I look at the layer in ArcGIS, there is no problem and all column labels and text is fine. Any idea what is causing the problem?


texts:
I have two polygons, a red one and a green one. These polygons overlap.

I would like to merge these two polygons. But when I merge them, I would like the attributes of the area where they overlap to be assigned the red polygon’s attributes, not the green polygon’s attributes. I do not want the area where they overlap to be assigned both polygon’s attributes.

i.e. When merging, I would like to assign a higher ranking to the red polygon.

I'm using ArcView 10.0


texts: What is correct link for adding GeoServer WFS to ArcMap?

texts:
What is the correct URL to enter in ArcCatalog WFS parameters when adding a WFS layer from GeoServer into ArcMap?

Here's a link to our GetCapabilities document:

https://weaveprod.ucdp.utah.edu/geoserver/ows?service=WFS&version=1.0.0&request=GetCapabilities

I am using ArcMap 10.1 and experiencing the same problem described here.  When I try to add WFS layers from GeoServer, "Feature Types" is blank.  That post says it is a bug with 10.1.  However, I added WFS from GeoServer not long ago, so I'm skeptical that the answer is to applicable to 10.2.  The only thing that changed was our GeoServer location.

I was able to bring in WFS layers from:

http://www2.dmsolutions.ca/cgi-bin/mswfs_gmap?version=1.0.0&request=GetCapabilities&service=WFS.

Although that is not GeoServer, I previously had GeoServer WFS working fine.

Here's an image of one URL I tried (which was:

https://weaveprod.ucdp.utah.edu/geoserver/sf/ows?service=WFS&version=1.0.0&request=GetFeature?

No matter how I change it, no layers are coming up in feature types.  There shouldn't be any problem with my WFS layers because I can preview them as WFS through GeoServer.


texts: Mapnik tiles have latitude distortion, projection issue?

texts:
I'm trying to work out a kink in setting up mapnik-node.

My geos are in postgis, and I allllmost have them rendering properly, but my GIS experience is lacking so I can't seem to determine why this isn't working correctly.

Here's an example overlayed on a map. I'm getting a weird distortion on the latitude but not on the longitude. http://cl.ly/image/0K3i1x3n2L3v

I can provide a snippet of the mapnik code, if that helps. But hoping this is just a simple bug—some kind of projection issue that I'm not seeing.

My polys are stored in the db with srid=4326.
The bounding box I calculate for the tile x/y/z matches perfectly with what should work.


texts: QGIS 2.2 Print Composer - pointer will not select or move items

texts:
I refer to the questions at QGIS 2.0 - can't select items in Print Composer and QGIS 2.2 Print Composer does not display map with some datasets

My question may be a duplicate, but the answer given in the first of those questions will not resolve the issue I am experiencing. If the CRS thing from the second one is a more likely cause, I would like a clear explanation.

I am using Windows 8.0 and the CRS of my map is EPSG 27700

It happens most often when I select my label, but it can be after selecting other items. Everything freezes on the layout - the selection pointer has no effect. I can add new items, but I can't move them or select anything else. I can select all or deselect all using the menu command, but still the pointer has no effect when clicked on the layout. I can lock and unlock using toolbar buttons, but not with a right click. I can use toolbar buttons to send selected items to back or front, but it makes no difference to the lack of pointer click functionality.

If I open a new print composer from my template, or if I close and re-open QGIS, the composer may work ok for a while, but then the same thing happens again. It was more intermittent, but today it is happening almost immediately every time.

I have disabled all my plug-ins, but it made no difference.


texts: How to get the "queried" value from attribute table and put it in a textbox?

texts:
I am new in VBA, I have a problem about querying an attribute table....:

I have a table of a layer (wish.shp) with the fields: owner_ID, Value, Area.

I made 3 textboxes in a userform, I input the owner_ID, in one textbox, I want to get the corresponding value and area and display in other textboxes.


texts: How to recreate the 'standard deviation' option for a raster layer in Mapserver?

texts: In mapserver is there any way to recreate the 'standard deviation' option available in Arc's symbology for a stretched raster layer?


texts: How do I split a shapefile efficiently into mutually exclusive sets?

texts:
I have a shapefile containing about 240,000 features. I want to split this into five files with no overlaps. My current approach is to write a non-trivial expression that defines the first subfile, save the selection, then delete the selected features. I then repeat this process. The problem is that deleting the records is taking 10-15 minutes for each step. In principle I could skip the deletion at each step, and then select the final subset as the negation of the OR of the earlier queries, but the likelihood of writing that query correctly seems quite low.

I suppose I could create a new field corresponding to each subset and set the value to 1 or 0 by expression, then set the last field as the negation of the OR of the fields. But I'm nearly done, so I'm not sure I'll be trying this.

For future reference, is there a better way to do this?