texts: How to clip very large shapefile in QGIS using ogr2ogr?

texts: I need to clip a large shapefile (600,000 KB) to another smaller shapefile (10 KB).  The vector -> geoprocessing -> clip command is taking very long.  Answers to similar questions say to use "ogr2ogr" but I am not sure where to start in installing ogr2ogr.  I have no Python experience either, not sure if that makes a difference.


texts: Matching csv of postal (zip) codes to another shapefile

texts:
I have a csv of survey results and each participant supplied their postal code. I also have a shapefile (polygon) that is the area of each postal code. I'm trying to get the csv answers (one point for every participant) to be displayed within the boundaries of the polygon.

I'm not concerned where inside the appropriate polygon just that they are within the bounds. I am unable to use geocoding at this moment and thought some sort of join/merge would work.


texts: How to reassign small polygons to its nearest neighbor in R (similar to gdal_sieve.py function)?

texts:
I created a spatial polygon file from a classified gridded dataset using the grid2polygons function in R. I would now like to take that polygon file and merge any polygons containing less than 6 pixels with its nearest largest neighbor. In the image link below, you can see an example of small polygons circled in blue. I would like those merged with the larger surround polygon. Alternatively, if there is a way to create a minimum threshold of 6 clumped pixels prior to the conversion to polygon stage, any suggestions for that method would be great as well.

There is a function in Gdal: gdal_sieve.py that can do exactly this (http://www.gdal.org/gdal_sieve.html) if you would like to get an idea of exactly what I would like to do. Does anyone know of a way to do this in R?


texts: format of raster calculator expression in ArcPy for ArcGIS 10.3

texts:
I'm trying to write a raster value inverter based on a raster calculator expression.

Article ID: 000006694 Link

The pure math would be:
(([elevgrd] - Z_Max) * -1) + Z_Min

Anyhow, I can't figure out how to write the expression. I started with .format and that failed with unicode errors, I then thought if I can get it out as a text string that might help, I copied the sting format (kinda) by exporting a model builder model of the populated geoprocessing tool..

But after 3 hours of swear at one line, I'm out of ideas..

    import arcpy
    from arcpy.sa import *
    InputRaster = arcpy.GetParameterAsText(0)
    OutputRaster = arcpy.GetParameterAsText(1)
    desc = arcpy.Describe(InputRaster)
    RasterPath = (desc.catalogPath)

    RasterPropsMINIMUM = arcpy.GetRasterProperties_management(InputRaster,"MINIMUM")
    MIN = (RasterPropsMINIMUM.getOutput(0))
    RasterPropsMAXIMUM = arcpy.GetRasterProperties_management(InputRaster,"MAXIMUM")
    MAX = (RasterPropsMAXIMUM.getOutput(0))

    expression=("\""+"\""+"\""+"("+"\"" +str(RasterPath) +"\"" +" - "+str(MAX) +" )*-1) +"+str(MIN)+"\""+"\""+"\"")
    arcpy.gp.RasterCalculator_sa((expression), OutputRaster)


texts: Inserting ground control points in ASP NASA?

texts:
I am using PLEIADES stereo images and  the Nasa software : Ames Stereo Pipeline (ASP) in order to generate a Data Elevation Model (DEM).

In addition to that, I have some ground control points and I would like to use them to improve DEM quality. How do I do it ?


texts: Seeking up to date Italian administrative boundaries?

texts:
I am after Italian GIS data - administrative boundaries such as:

Regions (Regioni)
Provincies (Provincie)
Municipalities (Comuni)
Districts (Quartieri)
City quarters (Contrade o Rioni)
Zones (Zone)
Suburbs (Frazioni)

and Postcodes.

I can find data some data to download - for example
http://www.dati.gov.it/catalog/dataset

http://www.igmi.org/download.php

But I want to make sure I have the most up to date data.

Who is the official supplier of this data?
Which government department manages this GIS data?
When the data is updated how will I know?

In Australia  there is a government department that manages this data and every 3 months we are are sent updates if we pay some $.
I am not sure if there is an equivalent in Italy


texts:
I am new to GIS, I am using PostGIS(2.3). I have two problems I am not clear with

Problem-1: I have list of locations(longitude, latitude, 4326), what is best/correct projection to use while clustering if I want to use unit in meters? I have looked around, people used 2950, with this when I project points back on map markers are going somewhere else, here is the query I am using

select
    array_to_json(
        array_agg(
            row_to_json(t)
        )
    )
from
    (
        select
            to_char(
                row_number() over(),
                '1'
            ) as id,
            ST_NumGeometries(gc) as count,
            ST_AsGeoJSON(
                ST_Transform(
                    gc,
                    4326
                )
            ) as points,
            ST_AsGeoJSON(
                ST_Transform(
                    ST_Centroid(gc),
                    4326
                )
            ) as center,
            ST_AsGeoJSON(
                ST_Transform(
                    ST_ConvexHull(gc),
                    4326
                )
            ) as border,
            avg_cost as "averageCost",
            min_cost as "minCost",
            max_cost as "maxCost"
        from
            (
                select
                    unnest(
                        ST_ClusterWithin(
                            ST_Transform(
                                geom,
                                2950
                            ),
                            500
                        )
                    ) gc,
                    avg( cost ) avg_cost,
                    min( cost ) min_cost,
                    max( cost ) max_cost
                from
                    user_locations
                where
                    region = 'sanjose'
            ) f
    ) t

Problem-2: I want to control clustering with respect to map zoom, so that I can change distance and minimum number of points in cluster, for this want to move to DBScan algorithm, it is not clear how to construct convexhull along with dbscan

Please help me


texts: Polygon to Raster: Cells with multiple values

texts:
I'm trying to use R to rasterize a number of shapefiles and I'm having trouble dealing with cells on the edge of shapes, in other words the if a cell is partially covered by an NA value and partially covered by a value associated with the shape. If any part of the cell is covered by a shape I would like the final raster product to have that value. For example in the image below I show the original shapefile (in green colors) and the result of a rasterize operation using R (in grey-scale). All of the red circles show areas where I would like the raster to have a value because some portion of the cell is covered by a shape.

I've tried using various functions with rasterize (max, mean, na.rm = T) and haven't had any luck producing the result I'm looking for. Below are links to the datafiles and my basic procedure.

raster: ftp://ftp.cdc.noaa.gov/Datasets/cpc_us_precip/precip.V1.0.mon.ltm.nc
shapeFile: http://droughtmonitor.unl.edu/data/shapefiles_m/USDM_20110802_M.zip

R is my tool of choice but I'm happy with solutions in Arc or QGis as well. With the caveat that any solution needs to be relatively fast (I have 900 of these and I would ideally like to get this done in less than a day) and batchable (I'm familiar with model builder in Arc)

##Load libraries
library("rgdal") #Handling Shapefiles
library("raster") #Handling Raster files
library("ncdf4") #For reading raster from .nc file

##Load in shapefile as spatialpolygon
myShape <- readOGR(dsn = "USDM_20110802.shp", layer = "USDM_20110802")

##Load in Raster
#This will throw a waring about the variable name it can be ignored
#The netCDF file is multilayer but we only need on for a mask
myRaster <- raster("precip.V1.0.mon.ltm.nc")

##Rotate Raster-orignial file use longitude of 0 to 360 we want -180 to 180 to work with the shapefile
#This will also throw an error which can be ignored
myRaster <- rotate(myRaster)

##Rasterize the shape file
#I've tried all the following
shapeRaster <- rasterize(myShape, myRaster, field = "DM")
shapeRaster <- rasterize(myShape, myRaster, field = "DM", na.rm = T)
shapeRaster <- rasterize(myShape, myRaster, field = "DM", fun = max)
shapeRaster <- rasterize(myShape, myRaster, field = "DM", fun = max, na.rm = T)
shapeRaster <- rasterize(myShape, myRaster, field = "DM", fun = mean)

##My goal is for this cell to have a value of 0 instead of NA
#in the image it's the top circled cell on the lefthand side
shapeRaster[24,215]


texts: What is the best way to do reverse geocoding with PostGIS?

texts:
Is there a way to do reverse geocoding with PostGIS? There's a lot of information on the Web but I found it very confusing...

I have imported osm data into my own Postgres with PostGIS database. Now I want to make geocoding and reverse geocoding with my database. I don't know if there are some Stored Procedures, or a project somewhere that takes care of my issue.

The next step that I want to do is make a webservice that can handle this kind of requests, but that is another question.


texts: Editing QGIS Composer Locked Layer Shapefiles

texts:
Is there a way to edit a shapefile which is locked as a layer in the layer and style within the composer?

For example I had a shapefile consisting of 50 ponds which were added to the composer atlas as the coverage layer. Each map within the composer might contain trees around ponds then another with location of fishing points, all independently locked and set as locked layers and locked styles.

If I then wanted to add another 50 ponds to the shapefile for ponds, keeping the same name, will the composer update when re-opened to show all 100 atlases’ without having to edit all the maps individually?


texts: Displaying the results of a multi-table query

texts:
I am trying to display the results of the following:

cartodb.createLayer(map, {
user_name: "tobagoborn",
type: 'cartodb',
sublayers: [{
    sql: "(SELECT L.the_geom_webmercator, L.landing_sites, F.first_name, F.last_name,C.number_of_species, C.cartodb_id FROM catch C, lsites L, fishermen F WHERE C.have_catch = 1 AND C.site_code = L.site_code AND C.fisherman_id = F.fisherman_id;}", // AND C.catch_date > now()::date";
    cartocss: '#catch {marker-fill: #FF0000;}',
    interactivity: 'landing_sites,first_name,last_name,number_of_species',
    }]
})
.addTo(map)
.on('done', function(layer) {
    var sublayer0 = layer.getSubLayer(0);
    cdb.vis.Vis.addInfowindow(map, sublayer0,
       ['landing_sites','first_name','last_name','number_of_species'],
       {infowindowTemplate: $('#infowindow_template').html()});
})
.on('error', function(err) {
    console.log(err);
   //log the error
});

I understand how to display the results of a single table query.  but need some assistance with doing the same with the SQL statement above.


texts: Loss of Quality in UAV Imagery after building pyramids in ArcMap?

texts:
I am using arcgis 10.1 to plot and print maps with a background with UAV imagery. When I insert an image (Tiff format), it loses quality.

The program prompts me to create pyramids, once this is complete, the image loses quality. What causes this? Is there a parameter which can be used to retain the original image quality?


texts: Setting QgsLayoutItemLabel position by map coordinates in QGIS

texts:
I am making a layout export function and I would like to create a QgsLayoutItemLabel at a given latitude and longitude. But in QgsLayout the position is relative xy (150,100).

Is there a way to convert map coordinates to layout positions?


texts: How can I control what is passed to the OnCreate() method when a command or tool is initialized?

texts:
In an ArcEngine application can I control what is passed to the OnCreate() method when a command or tool is initialized?

(I am using a synchronization object because I have both a map and a pagelayout in my form.  Rather than just passing the active control I would prefer to just pass the synchronizer because it contains all the controls.)


texts: Adding graphic to graphicsLayer on top of feature layer using ArcGIS API for JavaScript 3?

texts:
I have multiple featureLayers and one graphicLayer for adding graphics. However I cannot add a graphic to this layer. Can someone help with what I'm doing wrong ?

featureLayer.on("click", function (evt) {
  attrs = evt.graphic.attributes;
  if(clickmode == 1) {
    var highlightSymbol = new esri.symbol.SimpleMarkerSymbol().setColor(new dojo.Color([255, 0, 0]));
    var location = new Point(evt.mapPoint.x, evt.mapPoint.y, evt.spatialReference);
    graphicLayer.add(new Graphic(location, highlightSymbol));
  }
}


texts: How to query with stored PostGIS time series raster data

texts:
I have time series raster data in PostGIS and I can't find a way to get the info as a time series from an XY point.

For example:

I would like to get the vale of the XY point in raster for the 13/08/2015 - 20/08/2015

I have one raster for each day in the year 2015, each raster has values of some data from remote sensing satellite and I want to see how the value is changing week on week.

and so on.


texts: For loop to convert .lyr files to KML files gives ERROR 000814?

texts:
I must convert .lyr files to Kml files. After searching, I already tested 3 codes.
Unfortunately, when running the 3 codes Python displays same error:

  Runtime error  Traceback (most recent call last):   File "", line 4, in    File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\conversion.py", line 2322, in LayerToKML     raise e ExecuteError: ERROR 000814: Invalid file type

I have no idea how to fix the error.

#1
arcpy.env.workspace = r"D:/ss3/layers"
ou2= r"D:\ss3\kml"
for fc in arcpy.ListFiles("*.lyr"):
    arcpy.LayerToKML_conversion(fc, 'buf_layer')
    arcpy.SaveToLayerFile_management('buf_layer', ou2 + os.path.splitext(fc)[0] + ".kml", "ABSOLUTE")

#2
arcpy.env.workspace = r"D:/ss3/layers"
ou2= r"D:\ss3\kml"
for fc in arcpy.ListFiles("*.lyr"):
    arcpy.MakeFeatureLayer_management(fc, 'buf_layer')
    arcpy.LayerToKML_conversion('buf_layer', ou2, "12000", "false", "DEFAULT", "1024", "96", "CLAMPED_TO_GROUND")

#3
import arcpy, os
arcpy.env.workspace = r"D:/ss3/polygons"
ou2= r"D:/ss3/layers/"
for fc in arcpy.ListFiles("*.shp"):
    arcpy.MakeFeatureLayer_management(fc, 'buf_layer')
    arcpy.SaveToLayerFile_management('buf_layer', ou2 + os.path.splitext(fc)[0] + ".lyr", "ABSOLUTE")


texts: QGIS Training manual 21.14.1 : Where is Query Builder in QGIS 3.6?

texts:
I am new to GIS and QGIS and I am going through the training manual. It says to "Open theQuery Builder by right-clicking on the all_terrain layer in the Layers panel, and selecting the Properties→Sourcetab."

This does not exist as far as I can tell in 3.6. No other posted questions have an answer. Help?

@ahmadhanb stated that I should do what is written out above. The problem I am facing is that this is what I get when I do as ahmadhanb suggests:

I am using windows 10 if that helps any.


texts: Back to QGIS 3.10 from QGIS 3.12 in Ubuntu 19.10

texts:
I've followed this steps to install QGIS in linux. In the past, it worked fine to install QGIS 3.10, which I want. Now, this code install QGIS 3.12. What should I change to get 3.10?.

# edit sources.list and add deb lines
sudo nano /etc/apt/sources.list

deb     https://qgis.org/ubuntu eoan main
deb-src https://qgis.org/ubuntu eoan main

#public key
sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-key 51F523511C7028C3

#install
sudo apt-get update && sudo apt-get install qgis qgis-plugin-grass


texts: Reprojection not working with raster even though properties say its in the correct projection

texts:
I have a raster file that keeps placing itself in the South Pacific Ocean, West of Peru. But it should be placing itself over Colorado. I have tried using the tool "Define Projection" in ArcGIS Pro and "Project Raster", but it still places itself in the ocean. I doubled checked the raster properties, and under Geographic Coordinate System it says NAD_1983_UTM_Zone_13N, but again it places itself in the ocean and not over Colorado.
Has anyone run into this issue before?
Does this information help? It's from the original file I have before I re-projected it. It should cover the entirety of Colorado.


texts: No SIP module on Mac OS after installing Boundless Desktop 1.1.0 product

texts:
After installing a Boundless Desktop 1.1.0 product, when I tried to start QGIS, it gave an following error and processing toolbox does not started.

  Couldn't load SIP module. Python support will be disabled.

  Traceback (most recent call last):   File "", line 1, in  ImportError:
  No module named sip

I tried to install SIP module using the pip in boundless command prompt was unsuccessful. I also taken care as given as answer. But that does not work. Please suggest the issues.


texts:
I made a little plugin in QGIS, and I got some strange results.
The program now mirrors the "Layers Panel". When I add or delete layers from the Panel the QListWidget module get the names of layers.
I connected the function to the legendLayersAdded and to layersRemoved SIGNAL.

THE PROBLEM:
The add part is ok, and if I remove ONE layer that's fine too. But if I remove multiple layers, the widget gets back the rest of layers' names multiple times (count of removed layers).

Community edit: For however many multiple layers are removed, the remaining layers appear in the list that many times. Ergo, if 3 layers are removed, and there is 1 remaining, that layer appears 3 times in the additional layer pane:

It's like the signal calls the function recursively!?
I tried with the layerWasRemoved Signal and the result was the same
I made a video of it because it is hard to explain.

[https://www.youtube.com/watch?v=NuLtANtlLtE&feature=youtu.be]

Here is the code:

def ablak(self):
layers = []
self.dockwidget.telepules_kezelo.clear()
layers = self.iface.legendInterface().layers()
print "Befor" # Befor layer list
print layers # Layer list "from Layers Panel"
print "After" # After layer list
i = 0 # Just a counter for debubing
lista = self.dockwidget.telepules_kezelo
for layer in layers:
    print layer.name()
    lista.addItem(layer.name())
    i = i+1
    print i # Counter for debuging
layers = []

QgsMapLayerRegistry.instance().legendLayersAdded.connect(self.ablak)

QgsMapLayerRegistry.instance().layersRemoved.connect(self.ablak)

I fade out the problem with transfer the self.dockwidget.telepules_kezelo.clear() line after the Print After line. And it's working fine, but I am really frustrated with it.


texts: How to build cached tiled services in ArcGIS API JavaScript

texts:
I'm looking for an example or guidance on how to build cached tiled services in ArcGIS API JavaScript. Basically, using layer {z}, row {y}, and column {x} data with my arcgis server. Like this:

https://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{z}/{y}/{x}

Which results in this: https://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/19/220428/141938

I can accomplish this using either OpenLayers or esri with Leaflet.js (samples shown below). I am leaning towards leaflet.js, and I may want to implement both for more advanced work later on, but for now I want to learn how to do this strictly in ArcGIS.

OpenLayers sample

<!DOCTYPE html>

    <title>XYZ Esri</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.6.5/css/ol.css" type="text/css">

  https://openlayers.org/en/latest/examples/xyz-esri.html

Leaflet.js with Esri  sample

<!DOCTYPE html>

    <title>Simple Leaflet Map</title>
    <meta charset="utf-8" />
    <link
        rel="stylesheet"
        href="http://cdn.leafletjs.com/leaflet-0.7/leaflet.css"
    />

  http://bl.ocks.org/d3noob/8663620


texts: sld parser for version 1.1.0 with Openlayers

texts: I am using Openlayers 2.11 and need to style a layer with a sld file version 1.1.0. the error it says is Uncaught Can't find a SLD parser for version 1.1.0. i have also tried it with Openlayers 2.12.


texts: Import > Feature Class (multiple) of ArcMap not loading all features

texts:
The Import > Feature Class (multiple) option of ArcMap is not loading all my features. When I do this it fails and only loads partial features (only 200k features instead of the full 800k). When I do Import Single Feature it runs with no problems. It is a basic shapefile import to a geodatabase feature class.
I have been doing this type of import for years with no problems and now out of nowhere it seems to be causing me trouble, slowing down my monthly data updates.
Why might this be happening on the import multiple features only?


texts: ArcPy expression syntax in CalculateField_management

texts:
I am trying create a new column in my shapefile by concatenating a fixed/defined value with the value in another column. I have a syntax error that I can't seem to catch. Below is my code, current output, and desired output. There is a similar question at Writing expression syntax for Calculate Field in ArcPy?, but my case is a little different.
I am using ArcMap 10.7 & Python 2.7.
myval = 1
arcpy.AddField_management(myshp, "newCol", "TEXT")
expression = str(myval) + '-0' + str("!otherCol!")
arcpy.CalculateField_management(myshp, "newCol", expression, "PYTHON_9.3")

Current Output

otherCol
newCol

2
-1

2
-1

2
-1

1
0

Desired Output

otherCol
newCol

2
1-02

2
1-02

2
1-02

1
1-01

I also tried expression = str(myval) + '-0' + str(int("!otherCol!")), but it threw the following error: ValueError: invalid literal for int() with base 10: '!feat_ds_id!'.
The field otherCol is type Double.

This syntax is closer, but still not quite there:
expression = '"' + str(myval) + '-0' + '!otherCol!' + '"'

otherCol
newCol

2
1-02.0

2
1-02.0

2
1-02.0

1
1-01.0


texts: How do I adjust transparency on heatmap so that I can see background raster map? QGIS 1.8.0

texts: I have created a raster heatmap for the first time, but I can't seem to adjust the transparency settings to let me see another background OS raster map beneath it. I have attempted to swap the positions of the maps around in the layers bar but this does not make a difference.


texts: 'infinity' PostgreSQL date/time input not showing in QGIS

texts:
I have a PostgreSQL table with a date-type column populated either with proper dates (e.g. 2016-01-01) or with PostgreSQL infinity date/time input.

When I display the table in QGIS 2.14.11-1 (via DB Manager), the rows containing infinity values display NULL values instead (see last row on screenshot below).

Therefore, I'm unable to apply a rule-based style for these rows that remain impossible to select via the following rules (both rules return 0 selected feature):

"date" IS NULL
"date" = NULL
"date" = 'infinity'

Any idea on this issue and how I could deal with it?


texts: Error with tool 'NoneType' object has no attribute 'boundingBox'

texts:
I have a bespoke tool that takes a point input with multiple house prices and takes an average for the prices contain within polygons. The tool usually works fine but we have used a different house price input which appears the same as what we have previously used. The error read:

Traceback (most recent call last):
File "C:/Users/*****/.qgis2/python/plugins\*****\dist_price.py",
line 51, in run
self.copy_features_in_catchment()
File "C:/Users/*****/.qgis2/python/plugins\*****\dist_price.py",
line 132, in copy_features_in_catchment
in_cmnt = self.check_pt_in_catchment(src_feat, 'A')
File "C:/Users/*****/.qgis2/python/plugins\*****\dist_price.py",
line 264, in check_pt_in_catchment
fids = index.intersects(pt_geom.boundingBox())
AttributeError: 'NoneType' object has no attribute 'boundingBox'


texts: Best choice for building static maps: PostGIS, SpatiaLite or Shapefile?

texts:
I usually build static topographic maps, which are mostly printed. I'm often using OSM and rasterfiles (mainly different DEMs). The areas I'm mapping seldom exceeds 20.000 km². I'm using exclusively open source software like QGIS, SAGA and GRASS.

I know the problems with shapefiles, but they are easy to handle. And Im not very familiar with postgis or SpatiaLite. I tried postgis, but got a lot of problems handling it, SpatiaLite seem to me more easy.

What do you think: Should I stay with shapefiles or switch over to databases like SpatiaLite or postgis. Whats the advantages and disadvantages in my case?


texts:
I'm trying to learn how to create Script Tools in ArcGIS Pro, so I thought I would start small with a tutorial on https://www.e-education.psu.edu/geog485/node/120

The tutorial has you make a simple buffer tool, but whenever I check the module it says invalid syntax.

Here is my script and the error message. I know that it is identifying "arcpy" on the last line as the problem. I just don't know why.


texts: ArcPy Changing Spatial Database Connection

texts:
We have moved some servers which had ArcGIS application and SDE on it to another server and the IP has changed. Now when I want to open my .mxd file it tries to connect to remote SDE with old connection settings and opening old works but really takes long time, consequently I want to change my database connection in .mxd file.

Is there any fast way to change SDE settings programmatically with not opening the .mxd file in ArcMap?

import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\Users\arc\Desktop\election2011.mxd")
brklist = arcpy.mapping.ListBrokenDataSources(mxd)
print brklist

Above code really takes along time. without this code i can open my .mxd file with same time above code and can change manually...


texts:
How can I efficiently subdivide one raster into multiple rasters using the arcpy.SearchCursor?

For example if I have a  slope and I was to subdivide into 3 rasters:

slope1 with values between 0 and 25degress

slope2 with values between 25 and 50degress

slope3 with values > 50degress

Here I should easily solve the problem with arcpy.AddFieldDelimiters

If for example i have one categorical file and I also want to do the same thing?
If the file I have has only one field with the identification of the category and it's a string.

Value    Count    Type
    1    1023     Protected Area 2000
    2    1000     Protected Area 2001
    3    9555     Protected Area 2002
    4    9500     Protected Area 2003
    5    200     Construction Site
    6    (...)    (...)

Should I use arcpy.AddFieldDelimiters?

Furthermore, if I want to aggregate certain categories in only one raster, other groups of categories in other rasters?

When using the Search Cursor we are doing the changes directly to the source file ?should I create a copy for each raster I want to produce?


texts: Loop through 2 folders and perform raster calculation on the same file name

texts:
I was able to run the below script from @EmilBrundage.

However, the raster to be calculated were not paired correctly.

I am trying to loop files from 2 folders and perform raster calculation afterwards. The sequence of calculation must be:

(1st) out-0225*out-0225 [from clip qx folder] + out-0225*out-0225 [from clip qy  folder]

(2nd) out-0226*out-0226 [from clip qx folder] + out-0226*out-0226 [from clip qy  folder]

(3rd) out-0227*out-0227 [from clip qx folder] + out-0227*out-0227 [from clip qy  folder]

and so on.

But then the script calculates like the image below.

The script used is the ff:

import arcpy, os, math
arcpy.CheckOutExtension("Spatial")

ws1 = r"E:\clip_qx"
ws2 = r"E:\clip_qy"
outws = r"E:\ras_calc_qxy"

arcpy.env.workspace = ws1

rasterlist  = arcpy.ListRasters()

for r in rasterlist:
    r1 = arcpy.sa.Raster(r)
    r2 = arcpy.sa.Raster(os.path.join(ws2,r1))
    result =  (r1*r1) + (r2*r2)
    result.save(os.path.join(outws, r))


texts: How to create polygon made from outside points of points shp?

texts:
I have DEM Data which is in SHP Format. After interpolating to raster, I need to clip the raster image. To do this, I need to transform the Point SHP into a Polygon SHP using only the outer-most points to create the polygon. A Convex Hull does not produce the required results as shown in the diagram - The red outline is the Convex Hull. The arrows show the areas which should also be clipped.

I Have FME and ArcGIS Advanced.


texts: Spatial Join works when invoked through ArcMap GUI, but "Out of Memory Error" through ArcPy

texts:
I want to join an input feature class (~11,000 polygons) with another feature class (~31,000,000 polygons). When I do this through ArcMaps 'Join - based on spatial location'  tool, it takes quite long but eventually succeeds (adding 'Sum' for all joined fields). I monitored memory consumption throughout and it seems 'healthy' to me (not even going close to my 64 GB RAM).

When I try to do the same with the follwing ArcPy code (as adapted from the documentation: https://pro.arcgis.com/en/pro-app/tool-reference/analysis/spatial-join.htm):

outfc = "VG250_GEM_WBUILD"
inputFeature = "VG250_GEM"
joinFeature = "buildings_all"

fieldmappings = arcpy.FieldMappings()

# Add all fields from inputs.
fieldmappings.addTable(inputFeature)
fieldmappings.addTable(joinFeature)

# set sum area field
areaFieldIndex = fieldmappings.findFieldMapIndex("Area")
fieldmap = fieldmappings.getFieldMap(areaFieldIndex)

# Get the output field's properties as a field object
field = fieldmap.outputField

# Rename the field and pass the updated field object back into the field map
field.name = "sum_area"
field.aliasName = "sum_area"
fieldmap.outputField = field

# Set the merge rule to mean and then replace the old fieldmap in the mappings object
# with the updated one
fieldmap.mergeRule = "sum"
fieldmappings.replaceFieldMap(areaFieldIndex, fieldmap)

# Delete fields that are no longer applicable
fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex("osm_id"))
fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex("name"))
fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex("code"))
fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex("type"))
fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex("fclass"))

arcpy.SpatialJoin_analysis(inputFeature, joinFeature, outfc, "#", "#", fieldmappings)

I obtain the following error message after around 10 minutes:

  Runtime error  Traceback (most recent call last):   File "",
  line 1, in    File "C:\tmp.py", line 259, in
  main()   File "C:\tmp.py", line 251, in main
  calc()   File
  "C:\tmp.py", line 105, in
  calc
  arcpy.SpatialJoin_analysis(inputFeature, joinFeature, outfc, "#", "#",
  fieldmappings)   File "c:\program files
  (x86)\arcgis\desktop10.7\arcpy\arcpy\analysis.py", line 533, in
  SpatialJoin     raise e ExecuteError: ERROR 000426: Out Of Memory
  ERROR 000426: Out Of Memory Failed to execute (SpatialJoin).

What's strange about it is that, again, as far I can tell from the task manager, it's nowhere close from actually running out of my memory.

So, I believe I'm doing exactly the same thing through ArcPy as through the GUI, yet getting no result. What might cause these differences? How to troubleshout this?

Edit: I'm using ArcGIS Desktop 10.7.1 and the following Python version is used:

  2.7.16 (v2.7.16:413a49145e, Mar  4 2019, 01:30:55) [MSC v.1500 32 bit (Intel)]


texts: How to remove polyline rectangle from the Canvas on the UncheckBox Condition in QGIS using Python?

texts:
def handleItemClicked(self, item):
        canvas = iface.mapCanvas()
        self.r_polyline = QgsRubberBand(canvas, False)
        #self.r_polyline = None
        if item.checkState() == QtCore.Qt.Checked:
            item = self.dlg.tableWidget.item(item.row() , 1)
            #QgsMessageLog.logMessage(str(item.text()))

            self._selectedChart.append(item.text())
            self.res = []
            respw = client.service.GetSQLQueryDataInJson(commandText="SELECT chart_minx , chart_miny , chart_maxx , chart_maxy from tblchartmaster where chart_id = '%s'" %item.text() ,  serviceUser="pmluser", servicePass="pmlservicepassdal")
            json_data = json.loads(respw)
            for items in json_data:
                MINX = items['chart_minx']
                MINY = items['chart_miny']
                MAXX = items['chart_maxx']
                MAXY = items['chart_maxy']
                #canvas = iface.mapCanvas()
                self.r_polyline = QgsRubberBand(canvas, False)  # False = not a polygon
                #canvas.scene().addItem(self.r_polyline)
                points = [QgsPointXY(MINX,MAXY), QgsPointXY(MINX,MINY), QgsPointXY(MAXX,MINY),QgsPointXY(MAXX,MAXY) , QgsPointXY(MINX ,MAXY)]
                #QgsMessageLog.logMessage(str([points].text()))
                self.r_polyline.setToGeometry(QgsGeometry.fromPolylineXY(points), None)
                self.r_polyline.setColor(QColor(255, 0, 0))
                self.r_polyline.setWidth(1)
                self.r_polyline.show()
                #canvas.scene().removeItem(self.r_polyline)

                #self.r_polyline.reset()
                canvas.refresh()
                self.finished = False
            for i in self._selectedChart:
                if i not in self.res:
                    self.res.append(i)
                   # canvas.scene().removeItem(self.r_polyline)

        elif item.checkState() !=QtCore.Qt.Checked:

            item = self.dlg.tableWidget.item(item.row() , 1)
            if str(item.text()) in self.res:
                QgsMessageLog.logMessage('pass')
                self.res.remove(str(item.text()))


texts: Merging two vector layers with nearby features in QGIS

texts:
I am using QGIS 2.18

I have two shapefiles containing bus stops

the first file(old_stops) has 144 stops, the second (new_stops) has 5.

There are 3 stops in the new_stops file that are practically the same as 3 of the ones in old_stops- so only 2 new stops have to be added. So if I combine them, I should get 146 stops.

However, being a point shapefile, and because data was collected by different users, the 3 stops in old_files lie just a few meters away from the corresponding one in new_stops, and also have slightly different names

Because of this, a union results into 149 stops.

Now because it is a question of just a few stops, I can delete stops manually. However, I want to know if there is an automated process to do this for posterity.

To summarise: I am looking for an operation that can combine two layers,whilst setting a tolerance level of a few meters. That way, if a pair of features from each layer lie within the tolerance, it will choose the features of one of the shapefiles(which can be specified)

The easiest way I can think of(for large data) is:
1.creating a buffer layer around the old stops,
2. then select by location of all the stops in new_stops that lie within/overlap/touch the buffer layer,
3 deleting those selected features in new_stops
4. union of the modified with the old_stops

Is there an easier way?


texts: Creating species richness map (as next stage after species distribution)?

texts:
I'm using ArcMap. I've read the messages related to this topic but they don't take me far enough.

Most of the advice is related to species distribution, but I need to create a species richness map. I need to be able to say, for example, that in grid cell P there are 5 species (a to e), and in grid cell Q there are 22 species, even though in P there may be 200 records for species a, 100 each for species b , c and d, and 10 for species e. I'm simply wanting presence/absence of species, then a count of presences in any one grid cell.

So:

I have created a fishnet 1 deg x 1 deg for Australia.
I have 5 shapefiles for the separate species, but I also have a single shapefile with all the species in it (AS IN THE IMAGE ATTACHED). Each individual record has lats & longs.

I have tried to merge the shapefiles and the fishnet by doing a spatial join with ONE_TO_MANY, but the resulting attrib table is not listing or counting concurrent species.


texts: Adding Text and defined placement to a label in SLD (from QML(QGIS))

texts:
I've created a beautiful rule-based Labeling format in QGIS 3.12.

1) I have added text to the labels.

for example using the value : 'Beton'
then the word Beton is added to the line

2) I also used placement, data defined, coordinate X, coordinate Y:

endpoint:
x(line_interpolate_point($geometry,$length-5))
y(line_interpolate_point($geometry,$length-5))
beginpoint:
x(line_interpolate_point($geometry,5))
y(line_interpolate_point($geometry,5))

When I save the layer style to SLD, it seems lot's of info is missing, and not being saved to the SLD.

In the picture you see the difference in the labels.

How can I add some 'own'text and placement to the label in SLD
(f.e. using the application Notepad++ and add text to it)?
Is there anyone, with experience, what to write/add in de Notepad-file?


texts: Printing a Form includes data grid and axmapcontrol

texts:
I have created a windows form with some controls in it like data grid and axmapcontrol, it seems that printform cannot print the content of the axmapcontrol. I need to create a table of data and a map referencing the table data, like printing map and attribute table in one page.

I'm looking for any suggestions.


texts: Working in Full Motion Video (FMV), finalized multiplexed video not present?

texts:
Working within Full Motion Video, I have a .mp4 video file and its corresponding .csv file (both from Contour 2+ camera). When video multiplexing, I am able to get a success message within ArcMap; supposedly resulting in a MISB compliant video.

But the finalized, multiplexed output video is not in the output location i selected into the Video Multiplexer! There is simply nothing there (no video file, metadata, etc.).

Does it save elsewhere, or is there a user error?


texts: How to add to the map all the images on a image collection with cloud pixels masked on Google Earth Engine?

texts:
I'm trying to add to the map automatically trough a function a set of images from a landsat8 Image Collection to which previously I applied some filtering and a cloud masking function. I'm able to add all the images using the approach commented here: Add/display all images of mycollection in google earth engine, but the problem is that this way the clouds are not being filtered. I can't understand how to pass the cloudless images to the addImage function created on the script, I think it has something to do with the server/client structure of the platform.
Here is the link to the code:
https://code.earthengine.google.com/?scriptPath=users%2Fpaullcfias%2Fmosaics%3Alandsat8_rev
And the code itself:
    /**
 * Function to mask clouds based on the pixel_qa band of Landsat 8 SR data.
 * @param {ee.Image} image input Landsat 8 SR image
 * @return {ee.Image} cloudmasked Landsat 8 image
 */
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate('2019-01-01', '2020-03-31')
                  .filterBounds(table)
                  .filter(ee.Filter.or(
                  ee.Filter.and(ee.Filter.eq('WRS_PATH', 8),
                  ee.Filter.eq('WRS_ROW', 60)),
                  ee.Filter.and(ee.Filter.eq('WRS_PATH', 8),
                  ee.Filter.eq('WRS_ROW', 61)),
                  ee.Filter.and(ee.Filter.eq('WRS_PATH', 9),
                  ee.Filter.eq('WRS_ROW', 60)),
                  ee.Filter.and(ee.Filter.eq('WRS_PATH', 9),
                  ee.Filter.eq('WRS_ROW', 61))))
                  .filter(ee.Filter.lt('CLOUD_COVER', 75))
                  .map(maskL8sr);

print(dataset);

var listOfImages = dataset.toList(dataset.size());
print(listOfImages);

var visParams = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};

Map.centerObject(table);

function addImage(image2) { // display each image in collection
  var image = ee.Image(image2.id);
  Map.addLayer(image.clip(table), visParams);
}

dataset.evaluate(function(dataset) {  // use map on client-side
  dataset.features.map(addImage);
});

Map.addLayer(table, {},"zona estudio");

var slider = ui.Slider();
slider.onSlide(function(value) {
  var int_value = value * (Map.layers().length() - 1) >> 0;
  Map.layers().get(int_value).setOpacity(1);
  for (var i = int_value + 1; i < Map.layers().length(); i++) {
    Map.layers().get(i).setOpacity(0);
}
});
print(slider);


texts: I have a KML file containing several hundreds of Polygon all around the globe, projected in lat/lon WGS84, decimal degrees. I need to enlarge all the polygons adding 23km at the top and at the bottom of each polygon. Is it possible to do it directly in QGIS without exporting everything in excel and work it out?


texts:
I'm trying to delete a layer/features in my WFS I don't quite understand how the GET request for this works. I've got this so far :

http://192.168.3.169:8080/geoserver/ows?SERVICE=WFS&VERSION=2.0.0&request=Transaction&operation=Delete&typename=NaturalEarth:ne_10m_bathymetry_A_10000

But this doesnt work. I simply want to delete a Layer from my WFS via a Get Request. Also how do I use this to delete/edit/update/insert features within a layer?

I've tried reading the docs but don't quite understand it , still new to WFS


texts: How to calculate average value within a polygon(in Qgis with Python)?

texts:
Now I have a point and polygon shapefile

and I want to calculate the points' average value within each polygons

How to do it with python?

I tried to find out but google gives me a lot of raster within polygon...

I just want to do it with shapefile


texts:
i am currently trying to build a geoprocessing application using OpenLayers and ArcObjects,
I have some features that user draws on openlayers and i need to send these features to arcobjects for further processing.so i need a common feature format between these two systems.

solution 1

first i have tried getting esri-json support on openlayers with agsjsadapter but i got some unhandled exceptions - still trying to succeed.

solution 2

currently i am trying to convert a geojson format feature to Igeometry or Ipolygon on arcobjects, but i couldn't find any support for this in arcobjects.

i guess i need to solve this issue with one of the solutions above.


texts: What is the size of a geom at a given scale?

texts:
I have a dataset with many overlapping grid cells. I would like to represent the cells in such a way that each grid cell is represented with one point. To avoid having the points overlapping each other I have chosen centroid fill. The centroid is a simple marker. I have figured that the best way of doing this is by adding a random offset to the centroid marker.

This will give a visualisation near to a hatched map as more overlapping grid cells will lead to more centroids filling the grid cell.
The expected format described in the expression string builder is:

string of doubles 'x,y' or array of doubles [x, y]

This expression works fine:
to_string((rand(-20,20)/10))||','||to_string((rand(-20,20)/10))

Seeing how it works I now want a dynamic solution. How can I replace the grid cell size in millimeters with a calculated per grid cell (geom). Is there a way where I can calculate the represented size on the fly based on a fixed scale or even better - dynamically when I look at the map in layout mode?


texts: Geodjango - User's uploaded shapefiles and model creation

texts:
I have a big interrogation about how it could be possible to have a GeoDjango app that can store in database the uploaded shapefiles (or other geospatial files) from users. When you load shapefiles in GeoDjango, you have to create a new model that matches the fields of the shapefile. How could it be possible to create a new model automatically each time a user uploads a files? You can't append the new model definition to your models.py each time? For me, that does not make any sense.

Does it means the best way is to use a generic model definition and only allow files upload that match this definition,and then to create a new instance of this generic model? Can someone explain me what is the best practice to create an app that can store users uploaded spatial files?


texts: ArcGIS Feature Layer's Definition query to change by scale

texts:
I have a feature layer of road polylines. I would like to select the type of roads displayed on my client as a result of scale level.

For example:

At scale 100000 I would like to return only roads from type "primary".
At scale 50000 I would like to return only roads from type "primary" and "secondary".
At scale 25000 I would like to return only roads from type "primary" and "secondary" and "tertiary".

I know how to write a query that does this - I only need to know how to reflect the scale in this query. Is there any way the SQL statement 'knows' how to read the current scale from the client?

I've been searching and found that it is possible through the use of Query Layer but I don't want a Query layer. I already have that one layer and that's all I need.


texts: Workflow for identifying gaps in polyline shapefile

texts: So I have a polyline shapefile showing the border of a city, however there are multiple gaps in the shapefile ranging from ~2m to well over 1km that are almost impossible to see when viewing the area even at a large scale.  I am trying to generate a polygon using this outline in ArcGIS, however due to the breaks I cannot.  Is there any way to at least identify which line segments are not joined?  I attempted this solution in python with no luck.


texts: NumPyArrayToTable() (SystemError: error return without exception set)

texts:
My goal for this script is to extract some information from an Overpass API query and ultimately join it with the table of another feature class I have created.

So far I have completed the following steps relevant to my problem:

Extract the desired data from the JSON dictionary
Put data into a pandas dataframe
convert said dataframe to a sorted numpy array, then back to nparray
export the data in my numpy array to an ArcGIS table via NumPyArrayToTable()

With an empty query the script seems to work and gives me an empty table as desired,
which leads me to believe, that at least syntactically my code is correct.

import requests
import arcpy
import pandas as pd
import numpy as np

output = r'C:\Users\...\Output\test'

# GET request to overpass API
osmurl = 'http://overpass-api.de/api/interpreter'
osmrequest = {....}

osm = requests.get(osmurl, params=osmrequest)

# Extract data from JSON dictionary nested in elements
data = osm.json()
data = data['elements']

# rearrange osmdata so only relevant attributes in 'tags'can be extracted as dataframe
for i in data
    ...
        del i['tags']

# Create a dataframe from sorted JSON dictionary containing desired data
df = pd.DataFrame(data, columns=['data_1', 'data_2'])

# coerce dtype 'O' object values to np array native float64 type
df['data2'] = pd.to_numeric(df['data2'], errors='coerce')

# remove  all rows containing NaNs (no data2 information)
df = df[~np.isnan(osmdataframe).any(axis=1)]

# Create a sorted np array
sorted_osmNParray = np.array(osmdataframe.to_records())

# reconvert to np array so it can be passed to NumpyArrayToTable method
osmNParray = np.array(sorted_osmNParray, np.dtype([('index', '<i8'), ('data_1', '<i8'), ('data_2', '<f8')]))

# Convert NP Array to table
arcpy.da.NumPyArrayToTable(osmNParray, output)

The second I plug in real data though the following error gets thrown:

SystemError: error return without exception set

The cause of this is the method NumPyArrayToTable().

My Numpyarray looks like this before passing it to said method:

array([(20536,   4054816, 6.), (20537,   4054817, 2.),
       (20538,   4239284, 3.), ..., (23731, 581630073, 5.),
       (23732, 587279616, 4.), (23733, 590872784, 7.)],
      dtype=(numpy.record, [(u'index', '<i8'), (u'data_1', '<i8'), (u'data_2', '<f8')]))

I have tried the latter dataframe to numpy array conversion with a variation (NaNs still in my data) and it results in the same error:

df = pd.DataFrame(data, columns=['data_1', 'data_2'])

df['data_2'] = pd.to_numeric(df['data_2'], errors='coerce')

# Convert df to np record array with df.values, then reconvert to unstructured np array
n = np.array(np.rec.fromrecords(df.values))

# make list of column names to replace default nparray ones
list_names = df.dtypes.index.tolist()
n.dtype.names = tuple(list_names)

# Create table
arcpy.da.NumPyArrayToTable(n, output)

Could someone please tell me what exactly I'm doing wrong here?
Do I have to maybe reference an already existing table?


texts: Create laspy file object without creating a new file

texts:
I would like to create a new laspy.file.File object in memory without having to create a reference .las file.

The basic procedure is

read sample.las
clip sample.las by converting the x and y coordinates into numpy arrays and processing, I now have a boolean mask of coordinates in an arbitrary polygon as True
apply this mask to the original set of coordinates and return a new laspy.file.File object

I would like to do this without having to open a laspy.file.File object in write mode, which seems to be the most recommended solution. Is there any way to achieve this?


texts: Extract water areas from aerial photogrammetry DEM in QGIS

texts:
We are working with peatland restoration sites providing aerial photogrammetry of the restoration progress. We want to differentiate between areas water and bare peat which can look very similar in the orthomosaic raster. However, in the DEM the water features are instantly recognisable as they have a jagged and irregular surface. We can hillshade the DEM to show these areas and manually map them, but is there an easy way to automate this in QGIS?

Started playing with OTB and object based classification but not sure how to achieve this from the DEM as it's a single band raster.

Update, photos show a section of the orthomosaic, then the classified raster with bare peat (black), water (blue), vegetation (heather - green), and unclassified (transparent, showing the ortho below). It would be a great help to make the water areas solid and export as a shapefile.


texts: How can I display world hillshade and a shapefile without losing original colors and hillshade details?

texts:
This is what I am trying to do, but I am using the World Hillshade file from ArcGIS online and a land status shapefile:

http://www.esri.com/news/arcwatch/0312/learn-a-new-method-for-displaying-hillshades-and-elevation-tints.html

I am fairly new to GIS, so I do not think I am asking google the correct questions. I basically want to overlay my land status shapefile over the world hillshade, but the only way I know how to display both is by changing the transparency on one or the other. This article seems like it will give me what I want.

I have started by creating a mosaic dataset. The only way I can think of proceeding is to somehow convert the world hillshade and the land status to raster so I can add it to the mosaic dataset, but I do not believe that is possible and now I am stuck... There is a lot I am missing to the conceptual side of GIS.

Is the world hillshade something that I want? What about a DEM, should I be using that instead? I do not really understand the difference between the two or if that even matters.... all I want is for my map to have nice mountainous features with the land status color on top, both unaltered.


texts: Using ArcGIS Online to calculate walking times from points?

texts:
I am currently planning an accessibility analysis for the recycling infrastructure in my city (waste collecting points) and I want to determine the walking time from each point so that I could seek for the areas which are not covered by the infrastructure.

Using a questionnaire, I've found out that the people in my research area usually walk up to 20 minutes in order to reach their collecting points, so I would like to build the "buffers" for 4 time intervals (less than 5 minutes, 5 to 10, 10 to 15 and 15 to 20 minutes).

What should I do in order to accomplish this task?

I seek any help (including tutorials, videos, reading materials and so on).

Below is an example of what I'm trying to get (the exemplified map shows the walking time from supermarkets)


texts: ModelBuilder models keep getting broken/corrupt

texts:
I've got a problem with my ArcMap 10.4 ModelBuilder:

Whenever ModelBuilder model is saved/exited/crashed/etc. it has the chance of getting corrupt in a quite weird way. The model can still be used as a tool, the only thing that is different is that all the tools that were incorporated in the model are changed from being yellow editable boxes to white lines, which can't be opened nor edited. The model becomes uneditable when this happens. It seemingly happens by pure random chance.

Sometimes this problem can be fixed by opening another model and opening the corrupt one again, or deactivating some extensions. The corrupted model then changes back to normal and can be accessed like normal. But sometimes it seems I can't get my models back to their normal layout.

Here is a picture of the effect of this bug:

The tools can't be opened/edited, while the results and input data can be edited.

Does anyone know a workaround to repair this problem/prevent this from happening?


texts: Writing PyQGIS script with operations on multiple set of rasters

texts:
I would like to perform a loop that takes from a given folder a set of rasters and for every raster makes the following operation, that in the specific case is the calculation of the monthly evapotranspiration, that consist in multiplying every raster of monthly average temperature for a given coefficient (that is specified in the two arrays C and Rg).
The code doesn't give any error but it doesn't work, because at the end of the loop appears what follow: Invalid layer, Gdal cannot open the dataset.
Does anyone have any suggestion?
import glob, qgis
import gdal
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer
C=[0.4, 0.37, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4,0.4]
Rg=[109, 158, 261, 376, 467, 506, 518, 448, 321, 207, 117, 94]
suffix_output=['PET_jan', 'PET_feb', 'PET_mar', 'PET_apr', 'PET_may', 'PET_jun', 'PET_jul', 'PET_aug', 'PET_sep', 'PET_oct', 'PET_nov' ,'PET_dec']
suffix_input=['Tjan', 'Tfeb', 'Tmar', 'Tapr', 'Tmay', 'Tjun', 'Tjul', 'Taug', 'Tsep', 'Toct', 'Tnov', 'Tdec']

outputpath = '/Users/macbook/Desktop/TESI/PROGETTO/'
inputpath = '/Users/macbook/Desktop/TESI/PROGETTO/raster T average/'

list=[0,1,2,3,4,5,6,7,8,9,10,11]
for i in list:
  inputrasterfile = QgsRasterLayer(inputpath + suffix_input[i] + ".tif")

  entries = []
  ras = QgsRasterCalculatorEntry()
  ras.ref = 'ras@1'
  ras.raster = inputrasterfile
  ras.bandNumber = 1
  entries.append( ras )

  new_path = outputpath + suffix_output[i]+ '.tif'

  evapotranspiration = str(C[i]) +'* ras@1 ' + '*(str(Rg[i]) + 50)' + '/(ras@1 + 15)'

  Pet = QgsRasterCalculator(evapotranspiration, new_path, 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries )

  Pet.processCalculation()
  iface.addRasterLayer(new_path)

Just to understand, why these statement doesn't work?
I tried to follow the former suggestion, but the software says that I can't multiply sequence by non-int of type 'str'. I just would like to understand the rule behind this calculation...
evapotranspiration = (str(k[i])  + '* ( 0.46  * ras@1 + 8')   * str(p[i]) + '-2)' *str(d[i])

Here I want to write this formula: [k(i)p(i)(0,46*raster+8)-2]*d(i)


texts: Identifying duplicate attributes in field using QGIS?

texts:
I have a point shapefile with thousands of points. It has an ID code field that is supposed to be unique. Every now and then the data entry clerk wrongly type the ID creating duplicates. Right now I'm manually scrolling the field to find the duplicate.

Is there another way to do this using the Search Query Builder?


texts: Issue using pyModis.createMosaicGDAL() on Debian

texts:
I'm currently trying to automate some data processing with MODIS data using pyModis but am running into an odd issue.
I've developed a script that I've run locally just fine. It downloads all the tiles for a date using downModis() and creates the mosaic using createMosaicGDAL() and then does some additional processing. Locally, the downloading takes around 30 minutes IIRC and the mosaicking takes around 45 minutes.
I've set up a Google Cloud VM to automate the processing with. When running on this machine, the downloading is much faster (10 minutes) but the mosaicking is taking incredibly long or is getting stuck.
My local machine is a MacBook Pro with 4 CPU cores and 16GB RAM and a 256GB SSD.
The VM is running Debian with 2 cores and 16GB RAM and a 100GB hard drive.
The image below shows the Disk I/O while processing. You can see a spike in red which is the reads and then it drops down. The writes then spike in blue but then plateau at 12 MiB/s for over two hours. The final mosaic is about 10GB so this can't be the writing of the mosaic (12 MiB/s * 2 hours = 86400 MiB).

I'm wondering if there's any chance the code could be getting stuck somehow in some sort of loop. I don't see why else it would keep writing data. Also, the disk space remains constant from the start at around 10GB so it's not just writing additional files or anything like that.
This might totally be an issue with GCP but I figured I'd ask here to see if I could get any help.


texts:
I have a geojson of multipolygons like so

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "ADMIN": "Aruba",
        "ISO_A3": "ABW"
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              -69.99693762899992,
              12.577582098000036
            ],
            [
              -69.93639075399994,
              12.53172435100005
...

Using Python, I would like to simply add more properties to every feature in the feature collection. I would like to add these properties from a Python dictionary.

So I would take dictionary values and write them to the properties in each feature.

I've been a little lost so far frankly.


texts: Why do two rasters with the same spatial reference not align?

texts:
I have two rasters which have the same spatial reference and cell resolution but they are slightly out of line?

Any idea why this is and how I can align them?

I am using ArcMap 10.4.1


texts: Removing columns in a SpatialPolygonsDataFrame in R?

texts:
My spatial polygon data frame (SPDF) contains too many columns (variables) and I want to remove most of the columns entirely.

I know how to do this with a regular data frame in R, but I am unsure how to do this when dealing with object of class SpatialPolygonsDataFrame?


texts: Can I publish a Bing map from OpenLayers Plugin in my thesis?

texts: I have QuantumGIS 1.8.0. I made a map that I'd like to use in my thesis. I made it using Bing road layer map. I found it in OpenLayers plugin. Do I have the rights to publish the map as a part of my thesis?


texts: Where can I download 64 bit geoprocessing for ArcGIS 10.3?

texts:
Where can I download the 64 bit geoprocessing for ArcGIS 10.3?

I know it can be done from the instalation CD, but I don't have access to it at my work.

Do I just have to install ArcGIS 10.1 Service Pack 1 for (Desktop, Engine, Server)
?


texts:
I'm trying to understand how using each part of an IP address influences the accuracy of geocoding. For example, if

15.201.xxx.xxx

is used instead of

15.201.123.xxx

what is the approximate distance in KM? I'm using the IP2Location Lite database to examine who has visited a website. The README states the following:

1. IP2Location Lite Edition is free package with accuracy up to Class C (192.168.1.X) only. It is restricted for non-commercial use.

However, I'm having difficulty finding documentation on what the accuracy of this actually means. Does this accuracy vary over time also?


texts: How to extract a flow chart from points connected with lines?

texts:
I have some points connected with arrows, like the following example:

How can I extract a flow chart like the following?

I tried to use the flowMapper plugin but it seems it is not what I want.


texts: Setting creation of feature dataset as precondition to iterative process using ArcGIS ModelBuilder?

texts:
When using ArcGIS ModelBuilder, how do you set the creation of a feature dataset as a precondition to an iterative process, and not have it re-create the feature dataset each time through?

This is the model that I'm using:

I need the model to first create the feature dataset Test_FD and then populate that with feature classes in TestGDB.  However, the model creates the feature dataset each time through the iteration.  How can I get it to only create the feature dataset once, and then run the iterative process?


texts: Empty network dataset creation in ArcGIS for Desktop?

texts:
I'm trying to create a new network dataset via ArcCatalog in a SDE database. I've imported a street featureclass into my dataset, go to create the network dataset with all the default options. After I hit go, it spins for a minute then comes back saying the creation is complete, but that "the new network dataset has been created no features were found within the extent of the network dataset."

What causes this and how do I resolve?


texts:
I was trying to update my carto.js code from v3.15.19 to v4.2. What I recognized that carto.js v4.2 is no longer supports Named Maps.
I researched a lot and checking into the documentation made me wonder that there is no Named Maps mentioned at all inside the documentation for V4.2.
So please can anyone confirm if there is still any way to add NamedMaps layers to the map using carto.js v4.2 or not?


texts:
I have downloaded Malaysia map from opensteetmap website and I have started to work on it in PostGIS. I need to find out the distance of some points to the nearest street over the map, so I used this query to find the distance in kilometer:

select name,ST_Distance(ST_Transform(ST_SetSRID(ST_Point(101.6638513,3.114944),4326),3857),way)*60*1.6 as d from planet_osm_line order by d asc limit 5;

name  | distance
"Lebuhraya Persekutuan";453.941853570964
"Lebuhraya Persekutuan";453.941853570964
"";1291.97891210748
"Lebuhraya Persekutuan";1741.0179885363
"Kerinchi Link";1990.39830038167

The results are strange and biased by the way. The street name is correct and they are the nearest street to the targeted point, however the calculated distance is wrong. I noticed that SRID in openstreetmap is 3857, so I converted my target point to that SRID. If you check these streets and this point over googlemap, you will notice that the distance between them is almost nothing. Do you have any idea why the calculation is so distorted? (I suspected that conversion between SRIDs has not done correctly)


texts: How can I make sure my layers are in the same CRS projection?

texts:
I have lots of layers involving output areas and buildings in a city, and a shaky grasp of what CRSs are.

Basically, I feel like I've clicked the default OSGB 1936 British National Grid at every opportunity it asked me to specify the projection, but it also keeps alerting me like "Input layers have non-matching CRS"

How would I make 100% sure that every individual layer is projected the same, and if they aren't, how do I fix that?


texts: Understanding why extracted contours offset in QGIS?

texts: I am a beginner at QGIS, and I am trying to extract 1m contours from 5m DEM data. I have clipped a small area from the DEM data (roughly 1.2km by 1.2km). When I go to Raster - Extraction - Contour, it gives the correct contours, but they are offset on the map by about 1km to the west, and about 140m to the south. Any suggestions on how to fix this?


texts: Calculate maximum and minimum elevation for a polygon layer from a DEM

texts:
I have a polygon layer and would like to know the min and max elevation within each polygon. My elevation layer is a digital elevation model. I'm not quite sure how to extract that information.
I am using arcMap 10


texts: gdal_calc raster calculation fails on 28 GB GeoTiff

texts:
I am trying to do a raster calculations on a 28 GB GeoTiff using gdal_calc.py (3.0.4) but the process fails
gdal_calc.py -A DEM_0.tif --outfile=DEM_1.tif --calc="A+0.33" --format=GTiff --type=Float32 --co="COMPRESS=LZW"

when the size of the output file becomes ~650 MB the command keeps running but the file stops increasing in size. It happens if I run the process on the SSD internally but also on an external SSD. Furthermore, I noticed that the rate of increase of the filesize is only 70MB / min (which is rather slow I think, although gdal_calc only uses one CPU).
Thanks to a comment below I also created another geotiff with tilesize 64x64. I ran gdal_calc.py overnight and pretty much the same thing happened - it was still running after more than 7 hours so I killed the process.
GDAL or other software solutions welcome.
shortened gdalinfo:
Size is 146140, 86076
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
NoData Value=-32767
Overviews: 73070x43038, 36535x21519, 18268x10760, 9134x5380, 4567x2690, 2284x1345, 1142x673, 571x337, 286x169, 143x85
Unit Type: metre


texts: Memory limit and computation time out for NDVI Map in Google Earth Engine

texts:
I am trying to get an NDVI time series over a specific area of the Everglades and have the below script. When I try to run it I keep getting a memory exceeded error and also a computation time out error. For the memory issue, I have tried to use the . Reducer function but it still does not work. Any help or suggestion are welcomed

https://code.earthengine.google.com/10a3ad7ac8ccaa10e5fbd9e6a841690f


texts: Saving values by a Python script in a PostgreSQL database

texts:
I try to save attributes in PostgreSQL/PostGIS
by a Python script.

So far I get a connection and I can set an additional attribute field
that I can fill with a value (it seems so).
But how can I save these values?

My script looks like:

from qgis.core import
QgsApplication,QgsVectorLayer,QgsVectorDataProvider,QgsDataSourceUri
from PyQt5.QtCore import *
from qgis.PyQt.QtSql import QSqlDatabase, QSqlQuery

uri2 = QgsDataSourceUri()
uri2.setConnection("server", "5432", "table", "user", "pw")
uri2.setDataSource ("public", "table", "geom")
vlayer2=QgsVectorLayer (uri2.uri() ,"ww_manhole","postgres")
print (vlayer2.isValid())
if (vlayer2.startEditing()):
    print ("ok startEditing")
else:
    print("not ok")

my_field_name2 = 'new_field'
if (vlayer2.addAttribute(QgsField(my_field_name2, QVariant.String))):
    print ("ok AddAttribute")
else:
    print ("not ok2")

my_field_value2 = 'Hello world!'
for f in vlayer2.getFeatures():
    print ("ok getFeature")
    f[my_field_name2] = my_field_value2
    if (vlayer2.updateFeature(f)):
        print ("ok updatefeature")
    else:
        print ("not ok updatefeature")

if (vlayer2.updateFields()):
        print ("ok updateFields")
   else:
        print ("not ok updateFields")
if (vlayer2.commitChanges()):
    print ("ok commitChanges")
else:
    print("not ok commitChanges")

The output shows that neither
updateFields works nor commitChanges.
What change should I make?

The output:

True
ok startEditing

ok AddAttribute

ok getFeature

ok updatefeature

not ok updateFields

not ok commitChanges


texts: how to draw a point(10,10) on QGIS raster layer?

texts:
I Have QGis complied code with mingw32 and i am using qgsmapcanvas class on QWidget, I am able to load the raster and vector layers properly on the widget i Qt.

Now, i want to draw some points on layer only not on map canvas, but whenever i try to draw points using qgsrubberband class let's say point (10,10), it draws onto the mapcanvas , but i want that point to be drawn on raster layer (10,10) position.


texts: Update a field from a PostgreSQL point table (wherever it's NULL) with certain field value from closest polygon?

texts:
Given a point table (table1) and a polygon table (table2) in PostgreSQL, I am trying to update a field from table1 wherever its value is NULL. This field should be updated with a certain field value (let's call it 'required_field') from the closest polygon in table2.

Both tables contain spatial index (gist) and are on projected coordinates.
There are about 1.8 million points and 11000 polygons.
97.5% of points have got field_1 = NULL.

To obtain the value of required_field from closest polygon:

SELECT t2.required_field
FROM table1 t1
CROSS JOIN LATERAL
    (SELECT required_field, geom
     FROM table2
     ORDER BY t1.geom <-> geom
     LIMIT 1) AS t2

To update field_1 from table1 wherever its value is NULL:

UPDATE table1
SET field_1 = ...
WHERE field_1 IS NULL

However, I don't know how I could combine these two queries to produce the result I need.

I am trying with the following query, but it seems to run forever (more than 4 hours and it keeps going)

UPDATE table1
SET field_1 = required_field
FROM table1 AS t1
CROSS JOIN LATERAL
    (SELECT required_field, geom
     FROM table2
     ORDER BY t1.geom <-> geom
     LIMIT 1) AS t2
WHERE field_1 IS NULL;

Is this query actually right? Any other way to do what I am after?

.......

explain throws the following:

"Update on midpoints_wal_scot  (cost=0.15..1318769591009.67 rows=15380552946848 width=1543)"
"  ->  Nested Loop  (cost=0.15..1318769591009.67 rows=15380552946848 width=1543)"
"        ->  Nested Loop  (cost=0.15..3429851.15 rows=3868432 width=2940)"
"              ->  Seq Scan on midpoints_wal_scot t1  (cost=0.00..81053.14 rows=3868432 width=38)"
"                    Filter: (sop_id IS NULL)"
"              ->  Subquery Scan on t2  (cost=0.15..0.86 rows=1 width=2902)"
"                    ->  Limit  (cost=0.15..0.85 rows=1 width=1447)"
"                          ->  Index Scan using sop_polyg_geom_idx on sop_polyg  (cost=0.15..7672.73 rows=11029 width=1447)"
"                                Order By: (geom <-> t1.geom)"
"        ->  Materialize  (cost=0.00..143643.71 rows=3975914 width=62)"
"              ->  Seq Scan on midpoints_wal_scot  (cost=0.00..81053.14 rows=3975914 width=62)"


texts: Postgis ST_AsGeoJSON() to openlayers security issues?

texts:
Im about to implement a python script serving my openlayers code with GeoJson data from PostGIS. However, here How to connect OpenLayers to PostGIS data?
I see " Unfortunately you cannot connect a web page directly to a database because of security concerns "

Could you please tell me about those concerns? Isnt using a proxy /server script good enough?

PS. I would ideally like to post in the above tread but lack the 50 reputation neccesary.. sorry about that.


texts: Automating conversion from .SHP to .GPX for thousands of shapefiles?

texts:
I am looking for a way to automate conversion from .Shp to .Gpx. I know how to save a shapefile into gpx in a layer, but it would take too much time to do 1 by 1 as we might receive thousands of shapefiles that need conversion.

I am not very prone in coding, especially in QGIS, so does anyone have an idea or maybe even a code to make this possible with QGIS?


texts: Is there a workaround to the version limitation for QGIS resource sharing?

texts: I am using Qgis version 3.0.2. After downloading the QGIS Resource sharing plugin (by akbargumbira) I noticed it's only supported for the versions 2.4-2.99. Is there a workaround to use it on 3.0.2 (haven't been able to find any comparable tools to import icons etc.), or am I better off downloading an older QGIS version?


texts: How to return elevation for point data in QGIS?

texts: I have a dataset for the Western USA and want to get elevations for the point data added to the point attribute table. What is the best way to batch process this so the whole operation is done at one pass?


texts: Shapefile list gone missing in dialog window

texts:
I have updated my script in ArcGIS 10.3 to include the copyFeatures_management geoprocessing tool because I want to only copy the features that already have the same projection, not reproject them;

#Loop through shapfiles in folder and reproject
 for fc in fcList:
     fcspatialRef = arcpy.Describe(fc).spatialReference.name
     if fcspatialRef != spatialRef:
         arcpy.Project_management(fc, outFolder + "\\" + fc, template)
     else:
         arcpy.CopyFeatures_management(fc, outFolder + "\\" + fc)

# Print shapefile Project results
     arcpy.AddMessage(fc)

But when I do, I loss what I have been trying to achieve with the shapefile list in the dialog window;

![enter image description here][1]

The script works fine without the copyFeatures_management included;

#Loop through shapfiles in folder and reproject
for fc in fcList:
    fcspatialRef = arcpy.Describe(fc).spatialReference.name
    if fcspatialRef != spatialRef:
        arcpy.Project_management(fc, outFolder + "\\" + fc, template)

# Print shapefile Project results
    arcpy.AddMessage(fc)

Any ideas please?

Indentation fixed, code working for me, but still printing both the projected and copied shapefiles names in list, when I only want the projected shapefile names printed.

@fluidmotion, the;

arcpy.AddMessage('{} was projected'.format(fc)) # add if you want projected fc printed to screen

isn't working at all, and,

arcpy.AddMessage('{} was copied only'.format(fc)) #add if you want copied fc printed

is printing both the projected and copied shapefile names in the loop.


texts: WMS or WFS for selecting features and show attributes?

texts:
I am beginner in MapServer and Openlayers.

I want to be able to select features and show their attributes by user click. I don't need to edit any feature.

I know this is possible with WFS, but don't know about WMS.

is this possible with WMS?

What are your suggestions?


texts: How to improve the cartographic design of graticules?

texts:
The attached image shows a conflict between the graticule and the state boundary on a map.  I have not found a way to display user-defined lines of latitude, which would solve the 42 degree line conflict.

How can I resolve this conflict in Arcmap 10.1 by having better control of where the graticule parallels are placed?


texts:
I am using OWSLib to retreive a single WMTS tile, and want to add a worldfile to that tile so I can present it in a more "simple tool". When opening both the calculated world file in Qgis the image was not positioned at the some point, which made me write everything down. Below can be used to solve the problem.
My input:
wmts.tilematrixsets['EPSG:28992'].tilematrix['00'].topleftcorner
(-285401.92, 903401.92)

wmts.tilematrixsets['EPSG:28992'].tilematrix['00'].scaledenominator
12288000.0

wmts.tilematrixsets['EPSG:28992'].tilematrix['00'].tileheight
256

wmts.tilematrixsets['EPSG:28992'].tilematrix['00'].tilewidth
256

wmts.tilematrixsets['EPSG:28992'].tilematrix['00'].matrixheight
1

wmts.tilematrixsets['EPSG:28992'].tilematrix['00'].matrixwidth
1


texts: OpenLayers: change cursor upon hovering over feature

texts:
I want to change the mouse cursor to a pointer when the mouse is hovering over a feature.  I currently have the following:

this.map.on('pointermove', evt => {
  if (!evt.dragging) {
    this.map.getTargetElement().style.cursor = this.map.hasFeatureAtPixel(this.map.getEventPixel(evt.originalEvent)) ? 'pointer' : '';
  }
});

It works, but I am wondering if there is a more efficient way of doing this such as having the browser take care of it via CSS.


texts: Using polygons to create heatmaps in QGIS?

texts:
I've been fiddling around with heatmaps on QGIS for a few days, but came across a small obstacle on my way. I have a dataset comprised of point data, corresponding to objects inside a 1 sq meter grid. Each of those objects has several data associated, namely their weight (in grams). So, I have already created a simple heatmap based on the number of points and another one using the grams value of each point as the weight parameter of the heatmap. This resulted in two different heatmaps, both of them interesting to my analysis.

Nevertheless, I have a different dataset that comes from the same grid - this dataset corresponds to objects that I know came from a specific square on the grid, but to which I don't have a specific point location (as I said, I only know they came from a specific square, but I don't have XY coordinates for them inside that square). I have all the other data for these objects, namely their weight and I wanted to use them to create a new heatmap, where I would attempt to represent the distribution of objects through the grid, based upon their weight. I was thinking of creating a new attribute on my grid shapefile, where I would insert the total object weight (both the ones to which I have XY coordinates and the ones to which I don't have that information) for each square.

I've been looking at the heatmap plugin on QGIS, but it only allows me to use point shapefiles as the input layer... Is there a way of creating a heatmap using polygons as the input layer? Remember, the polygons do not overlap, so I cannot use the approaches sugested elsewhere, regarding overlapping polygons.

I have already created a "points per polygon" shapefile and I guess I could also create a "total weight por polygon" shapefile, but I'd rather have a heatmap, since it would give out a less "boxy" visualization.


texts: PostGIS -> LineString + PointZ = LineStringZ

texts:
DATA

I have 2 tables in PostGIS:

linestable with MultiLineString geom
surveypointtable with Point + elevation attribute

GOAL

The goal is to rebuild the lines as 3D replacing line's vertices with the closets surveypoint.

LOGIC

Rebuild points as PointZ and get the collection of all points

ST_Collect(ST_MakePoint(ST_X(geom),ST_Y(geom),elevation_att))

Get line's vertices as points using

ST_DumpPoints(geom).geom

Find closest surveypoint using

ST_ClosestPoint(surveypoints.geom, linepoints.geom)

Check if the closestpoint is within 1m using

ST_DWithin(surveypoints.geom, linepoints.geom, 1)

Build line using only the closets surveypoints that are within 1m tolerance.

ST_MakeLine(ST_Collect(closestpoints.geom))

PROBLEM

When I use ST_ClosestPoint(p.3) it drops the z-index of the point and as a result I end up with a 2D line.

When I use ST_3DClosestPoint it takes wrong points - i.e. if I had point A 1m away in 2D from linepoint but at the same elevation and point B at the exact same coordinates as linepoint but with 10m difference in elevation it would take point A.

CODE

    DROP TABLE IF EXISTS newlines;

    CREATE TABLE newlines AS
    SELECT gid,
          ST_SetSRID(ST_MakeLine(geom),27700)
    FROM
      (SELECT gid,
              path,
              ST_AsText(
                ST_Collect(
                  closestsurveypoint)) AS geom
      FROM
        ( SELECT linepoints.gid AS gid,
                  linepoints.path AS path,
                  linepoints.geom AS linepoint,
                  (CASE ST_DWithin(
                        ST_SetSRID(surveypoints.geom, 27700),
                        ST_SetSRID(linepoints.geom, 27700),
                        1)
                      WHEN TRUE THEN ST_AsText(
                        ST_ClosestPoint(
                          ST_SetSRID(surveypoints.geom, 27700),
                          ST_SetSRID(linepoints.geom, 27700) ))
                  END) AS closestsurveypoint
          FROM
            ( SELECT gid,
                    color,
                    (ST_DumpPoints(geom)).path AS path,
                    ST_AsText((ST_DumpPoints(geom)).geom) AS geom
            FROM linestable) AS linepoints,

            ( SELECT ST_AsText(
              ST_Collect(
                ST_MakePoint(
                  ST_X(geom),
                  ST_Y(geom),
                  elevation))) AS geom
            FROM surveypointstable) AS surveypoints ) AS newlinepoints
      GROUP BY gid,
                path
      ORDER BY gid,
                path) AS newlines
    GROUP BY gid


texts: Merging several DEMs into one using ArcGIS Desktop?

texts: Is there any ArcGIS 10 tool to merge several DEM into a single one?


texts:
I am trying to generate a hot-spot in crime profiler based on the total amounts of cost associated with a particular location.

I have multiple locations but I do not want a standard hot-spot showing the spread of the locations/incidents. Instead I want to use the sum of costs associated at the locations to generate the hot-spot.

In my data I have

Location
X and Y
Sum of minutes
Sum of costs
Sum of incidents

How do you do this?


texts:
I have a point feature that I want the closest address for, I have done this by a spatial join on the closest known address point.

However, I have the address now for each point but because of the way the address point feature that I used for the spatial join was set out the address could appear in 1 of 2 columns (see screenshot)

Is there a simple way I can just bring all the addresses into 1 new column from the 2 address columns, so I have a list of addresses for each point?

At the moment the addresses are scattered between 2 columns which makes it very hard to work with.

(I marked out in red the addresses for privacy reasons)


texts: Create lines by id from two sets of points

texts:
I have two sets of points and I want to connect them to form lines by the id of the points. I.E. I have the following feature classes of points:

FC_1

ID EASTING NORTHING
A1 222222  222222
A2 333333  333333

FC_2

ID EASTING NORTHING
A1 222822  222222
A2 333933  333333

I want to create a line feature class that joins by ID, so one line for ID A1 and another line for ID A2. How can I achieve this?


texts:
I am assisting on a project in the Philippines to map some areas. I've been given a sheet the has points, but its unlike anything I've seen before and have no idea how to put it in a format that would allow me to upload it into QGIS.

A brief example of what I am working with is like this:

"Beginning at a point mark '1' on Protected Area Boundary Map, being N 49-11-05 E, 5,221.294 meters from BHL-No.3080 (PRS 92 Control point: Latitude 09-57-30.57, Longitude 123-57-19.69) of the Tubigon, Province Bohol

Thence          N 81-30-40 E     203.226     m. to corner     2 ;
 Thence          S 60-20-27 E     248.566     m. to corner     3 ;

Thence          S 47-13-20 E     237.059     m. to corner     4 ;

and so on and so forth it continues like this until returning to point 1.

Looking for any insight on how to put this in a format in which I can upload these points to QGIS.


texts: "Error in .local(.Object, ...) : Dataset copy failed" in r raster

texts:
I'm trying to write raster bricks using lapply but I get the following error:
Error in .local(.Object, ...) : Dataset copy failed
Calls: lapply ... <Anonymous> -> new -> initialize -> initialize -> .local
Execution halted

Here is my entire code:
    library(raster)
    library(dismo)
    library(parallel)

    setwd("/scratch/bntjoa002/chelsa_cmip5_ts/")

    # set the number of cores
    #detectCores()
    mc.cores=80

    #get model names
    #modelnames<-c("ACCESS1-3", "MIROC5", "CESM1-BGC", "CMCC-CM") #could also read a csv file of the model names
    modelnames<-c("ACCESS1-3")

    #experiment names
    experiments<-c("rcp45", "rcp85")

    #specify years
    start.year<-2070
    end.year<-2089

    # start the loop to calculate bioclim vars

    for (i in 1:length(modelnames)){

      model = modelnames[i]

      for(j in 1:length(experiments)){

        experiment<-experiments[j]

        tasmax<-brick(paste0(getwd(),'/tasmax','/CHELSAcmip5ts_tasmax_',model,"_",experiment,"_",start.year,"-",end.year, ".nc"),
                      varname='air_temperature')
        tasmin<-brick(paste0(getwd(),'/tasmin','/CHELSAcmip5ts_tasmin_',model,"_",experiment,"_",start.year,"-",end.year, ".nc"),
                      varname='air_temperature')
        pr<-brick(paste0(getwd(),'/pr','/CHELSAcmip5ts_pr_',model,"_",experiment,"_",start.year,"-",end.year, "_V1.1.nc"),
                  varname='precipitation_flux')

        # make a vector with the years in climate data time series
        years<-seq(start.year,end.year,1)

        nyears=length(years)

        # make index for the year
        ind=sort(rep(1:nyears,12))

        # Calculate bioclim variables for each year using mclapply
        bcbricks=mclapply(unique(ind),function(x){
          #extract just the climate data for each year from the raster brick
          tasmax.year=tasmax[[which(ind==x)]]
          tasmin.year=tasmin[[which(ind==x)]]
          precip.year=pr[[which(ind==x)]]
          # calculate bioclim variables
          biovars(precip.year,tasmin.year,tasmax.year)
        },mc.cores=mc.cores)

        #if you want to save all bioclim vars for a year as a stack for each year
        #For FUTURE climate projections, stack all bioclim variables for each year
        # if(experiment=="rcp45"| experiment=="rcp85"){

        # output directory
        #   outDir=paste0('/ClimateModelData_NASA_GDDP_Bioclim/future/',model)

        #   lapply(seq_along(years),function(x){
        #     writeRaster(bcbricks[[x]],file=paste0(getwd(),outDir,'/',"Bioclim1-19","_",model,"_",experiment,"_",years[x],'.tif'))
        #   })
        # }

        #if you want to save one bioclim var as a stack across all years
        if(experiment=="rcp45"| experiment=="rcp85"){

          # output directory
          outDir=paste0('/scratch/bntjoa002/chelsa_cmip5_ts/',model)

          # stack individual bioclim variables across years
          lapply(1:nlayers(bcbricks[[1]]),function(x){

            bcvar.by.allyears=stack(lapply(bcbricks,function(y){ y[[x]]}))

            writeRaster(bcvar.by.allyears,file=paste0(getwd(),outDir,'/',names(bcbricks[[1]])[x],'_',model,'_',experiment,'_',start.year,'_',end.year,'.tif'))
          })
        }
        print(j)
      }
      print(i)
    }

I'm submitting the job to a cluster, and there is definitely enough space on the drive. Any ideas what the issue might be?


texts: Can't load WMS service, "Failed to parse capabilities"

texts:
Using QGIS 3.8.1 and trying to load the QLG Gov base imagery WMS layer. It was working perfectly fine but then suddenly just stopped working and gave the message "Failed to parse capabilities". Tried changing the CRS's, but that didn't help.

The server I'm trying to load is "https://gisservices.information.qld.gov.au/arcgis/services/Imagery/QldBase_AllUsers/ImageServer/WMSServer"

Anyone know what's going on?

You can probably tell I'm new to using GIS software, but willing to learn and get there!


texts: ArcPy Update Cursor after nested Search Cursor giving tuple assignment error?

texts:
I am running ArcMap 10.4.1. I am running a script that uses a DA search cursor nested inside an update cursor to assign variables from values in a table. The script then writes those same variables to a feature class in the update cursor. I have done this many times in the past, but now I am getting a tuple assignment error when I try to assign row values in the update cursor. I know tuples aren't mutable and I am not trying to assign values in the search cursor. What is the problem here? This syntax has worked many times in the past.

The script fails on the first attempted row assignment after the search cursor:

    row[2] = tblAppAdd
    TypeError: 'tuple' object does not support item assignment

Here is the script portion (edited for length):

with arcpy.da.UpdateCursor(u'lyr_patientsPts', fcFields) as fieldsCursor:
    for row in fieldsCursor:

        appID = unicode(row[1])

        # This is a custom function I found to simplify writing SQL search statements
        # for the where_clause in the search cursor.
        tblWhere = WhereClauseFromList(tblPatients, u'APPLICATION_ID', [appID])

        # Get values of the fields in the table to put into the points feature class
        with arcpy.da.SearchCursor(tblPatients, tblFields, where_clause = tblWhere) as tblCursor:

            for row in tblCursor:
                tblAppAdd = unicode(row[2])
                tblAppCity = unicode(row[3])
                tblAppState = unicode(row[4])
                # ... etc there are ten more fields

        # Write values from table to points feature class
        row[2] = tblAppAdd
        row[3] = tblAppCity
        row[4] = tblAppState
        # etc

        fieldsCursor.updateRow(row)


texts: Problem with encoding UTF-8 chars in Windows with PyQGIS

texts:
I have the problem to show characters (äöü) in a utf8 format in Windows. In Linux there is no problem and everything is displayed correctly. I use at the second line:

# -*- coding: utf-8 -*-

and further for importing and opening the shapefile:

base = os.path.basename(str(filename))
self.vlayer = QgsVectorLayer(filename, base, "ogr")

layer = self.vlayer
provider = layer.dataProvider()
layer.select(provider.attributeIndexes())
addA=[]
for elem in layer:
    attrs_r = elem.attributeMap()
    for (k,attr) in attrs_r.iteritems():
        addA.append(attr.toString())

I also tried:

addA.append(u'%s'%attr.toString())

and

layer.setProviderEncoding(u'UTF-8')
layer.dataProvider().setEncoding(u'UTF-8')

But it has no effect and the characters are still not displayed correctly. Is there somehow a possibility to assign directly utf-8 to the shapefile?


texts:
I would like to get a list of modifications on the OSM database, in a given bounding box, between two given dates and impacting a given set of tags.

Does anyone know how to do this?

Any command-line tool, webapp such as Overpass Turbo, or query API such as XAPI?

I have seen it's possible to get history of a map by its location and zoom level (example: http://www.openstreetmap.org/history#map=16/45.9605/5.3391 ), but this return groups of modification that contains modifications outside of the current map).

I have also seen a history browser  to browse history of a given object (node, way or relation).


texts: Adding features to PostGIS layer using PyQGIS?

texts:
Somehow I can't manage to add a feature to my postgres layer. Is this only possible in edit mode?

My layer

>>> cp
<qgis._core.QgsVectorLayer object at 0x7f541c0e0348>

>>> cp.source()
u'dbname=\'work_viewtest\' host=myhost port=5432 user=\'b\' password=\'bxxxxxx\' key=\'id\' table="public"."v2_calculation_point" (the_geom) sql='

This is how I create my feature

f = QgsFeature()
f.setGeometry(QgsGeometry.fromPoint(QgsPoint(123, 456)))
f.setAttributes([10, 100, "ref", -1])
cp.dataProvider().addFeatures([f])
(True, [<qgis._core.QgsFeature object at 0x7f539ffa33e0>])
cp.updateExtents()
cp.updateFields()

The addFeatures() method returns True but I don't see any features in my layer, nor in the database table itself.


texts: Conversion huge raster (50 000 000 cells) in to .xyz

texts:
I do have DEM approx. 40x5km with cell size 2x2m, which gives approx. 50 000 000 points.

I need to export it to .xyz (1pt = 1 line in format x y z) file in order to use it in MIKE21.

It seemed like there is loads of tools that will do it (I do have access to Arcinfo 9.3.1 and 10), though I'm stuck for 2 days already. I tried some solutions and all have failed.

I did:
1. SA/EXTRACTION/SAMPLE, receiving table, saving it as .dbf, and trying to import it in to Access 2007. I could not import .dbf produced in Arc into Access. Exporting table to .txt gives me always error.

Raster Calculator and Sample(). This worked great. Raster converted in 5 minutes producing 1.2GB file, but in order: ZXY not XYZ. MIKE do not allow to change order of columns, so I need to change order of columns. No clue how to do it. Tried Access 2007 again, but failed on import data.
Script Raster2XYZ ver. 1.1 from ESRI website (in VB). Can not install it properly now (gives an error), I used it before and it worked 6 months ago. However I remember that it worked really slow back then. I am not much of a programmer and can not find the bug.
Grass GIS and r.out.xyz have bug. Problem with double or float writing(can't remember now).
I have R script which does export .asc to .xyz. It take 2.5h to export 40 000 000 cells. Output file is 4.3GB. MIKE will not swallow it.

Any ideas how to export such raster to .xyz? Or how to reorder file from pt.2 to receive XYZ instead ZXY?

I always used Access to convert data (ie. .dbf to .txt(.xyz)), but now it seems like no option. Is there any other solution?

I do have HP workstation Z600 (Xeon 2.4Ghz, 8GB RAM, GF Quadro FX 3800).

Any sugestions?


texts: How to make new raster with higher values than the older one?

texts: I'm working in ArcGIS and I have an elevation model of seaside. What I should do is to model an uplift of land, so make a new raster dataset in which will be higher values (for example to make a new dataset where all the values for land, not for sea, will be 1 meter higher). Please, do anybody know how to do that?


texts: Manually changing attribute values for PostGIS layer in QGIS?

texts:
I am trying to edit the attributes of a PostGIS layer. I am able to create a new field using field calculator, but when I try to change a specific cell's value following documentation:

  typing the new value directly in the cell, whether the attribute table
  is in table or form view. Changes are hence done cell by cell, feature
  by feature;

, nothing happens when I type in the cell.

I've enabled the edit mode (and also tried the multi-edit) mode. This is with QGIS 3.0.1.

Are manual edits allowed for PostGIS layers or am I missing something?


texts: Rearranging a list of coordinates within a list of pair simultaneously while running a php code

texts:
The question is related to rearranging a list of coordinates within list of pair while running a php code simultaneously. I am new with php language and managed somehow to parse a kml file in php. The file came out with the following output:

Endeavour Hills
145.2355671,-37.9742105 145.2353525,-37.9749209 145.2343225,-37.9789804
Narre Warren
145.2355671,-37.9742105 145.2353525,-37.9749209 145.2343225,-37.9789804
Cranbourne
145.2355671,-37.9742105 145.2353525,-37.9749209 145.2343225,-37.9789804

(There were more coordinates for each records and I tried to keep it concise)

So my code could parse the coordinates from kml file. But I wanted the coordinates get rearranged simultaneously during the code run. For example it should look like following.

Endeavour Hills
145.2355671,-37.9742105
145.2353525,-37.9749209
145.2343225,-37.9789804
Narre Warren
145.2355671,-37.9742105
145.2353525,-37.9749209
145.2343225,-37.9789804
Cranbourne
145.2355671,-37.9742105
145.2353525,-37.9749209
145.2343225,-37.9789804


texts: Creating table containing contents of directory with FME 2014?

texts:
I have to match the contents of a directory with the attributes of a FeatureClass. The FeatureClass field contains the names (without file endings) of documents which are either .pdf or .jpg files.

I have two directories - one for PDFs and the other for JPGs. I need to match the PDF file names with the FeatureClass Attribute and if they match rename them with the PDF file ending. The ones that are not matched are then obviously JPG.

How can I create the list of files using FME2014?


texts: Can't save any polygon layer: getting strange error message

texts:
After using QGIS 3.12 for some time, I ran across a problem where I could not save any polygon that had a text field. Even if I create a new project with nothing in it, make a polygon, and add a text field, I get the following error when I try and save the layer.

I tried deleting the QGIS 3.12 directory (I downloaded it with the web installer) and installing QGIS 3.10 through the standalone installer, but the same issue comes up even on a fresh install! Is there some kind of temporary directory used by the program that I need to clear as well?


texts:
How can I create mapfile (.map) for MapServer 7.0 in QGIS 2.18 or 3.0?

Or are there any tools which are part of the core MapServer distribution?

And what is the best visual tool for create mapfile now?


texts: Interpolate temperature data (shapefile) using R

texts:
The answer below solved my fist question, so I edit my question as follows:

I have a shapefile given with the attributes lat, lon, elevation and temperature.

So far I have a script that interpolates only the temperature data:

#read shapefile:
shape <- readOGR("C:/Users.../Testarea.gdb", "temperature_data")
#grid shapefile:
shape.grid <- spsample(shape, type = "regular", offset = (c1,1))
gridded(shape.grid) <- TRUE

#kriging:
test <- krige(Temperatur~1, shape, shape.grid)
spplot(test[[1]])

Now my problem is that I should also need to consider the elevation, since the temperature is often dependant on this parameter. How can I do this?


texts: Iterate polygon column in ArcGIS model builder

texts: I have a polygon with many columns. I want to convert into raster. I could do it one by one. But I want to iterate for each column of the polygon by model builder or any other way. Here you can see the attached table of the polygon. For example I want to make one raster for Abst85, another one for Abst91, .....Then save them in separate


texts: Field Calculator returning incorrect length and area values when QGIS project is reopened

texts:
Similar to Field Calculator returning incorrect length and area values, but as always slightly different.

I have QGIS 2.0 and when I create a polyline shapefile, the field calculator returns the correct lengths.  However the next time I open the project or the shapefile in a new workspace, when I click and the field calculator to update the line lengths, the numbers (in one instance) go to the 10,000,000's.

I have checked the project workspace it is in UTM Z51S and the units are in metres and the shapefile projection is the same.

I have updated the values to be zero first and then use $length to no avail.  I have even tried to changing the projection to WGS84 with units as metres and the values are in degrees.  It even happens when I create the values in a new field.


texts: I routinely use a two-monitor setup with my portable Precision workstation.  When I'm in laptop only mode, my Table of Contents may appear briefly as a shadow on the extreme margin of the screen, then it disappears and I am unable to view it or the catalog windows.  I'm using ArcMap 10.0 and Windows 7.0