texts: Conditional simulation (with Kriging) in R with parallelization?

texts:
I am using gstat package in R to generate sequential gaussian simulations. My pc have 4 cores and I tried to parallelize the krige() function using the parallel package following the script provided by Guzmán to answer the question How to achieve parallel Kriging in R to speed up the process?.

The resulting simulations are, however, wrong, it looks a geometry problem, but i can't find how to tackle the problem. Next i will provide an example (4 cores). You will see that after running the code, the maps from parallelization have vertical lines, and are different from the ones using only one core at the time.

library(gstat)
library(sp)
library(raster)

# create a regular grid
nx=100 # number of columns
ny=100 # number of rows
srgr <- expand.grid(1:ny, nx:1)
names(srgr) <- c('x','y')
gridded(srgr)<-~x+y

# generate a spatial process (unconditional simulation)
g<-gstat(formula=z~x+y, locations=~x+y, dummy=T, beta=15, model=vgm(psill=3, range=10, nugget=0,model='Exp'), nmax=20)
sim <- predict(g, newdata=srgr, nsim=1)
r<-raster(sim)

# generate sample data (Poisson process)
library(spatstat)
int<-0.02
rpp<-rpoispp(int,win=owin(c(0,nx),c(0,ny)))
df<-as.data.frame(rpp)
coordinates(df)<-~x+y

# assign raster values to sample data
dfpp <-raster::extract(r,df,df=TRUE)
smp<-cbind(coordinates(df),dfpp)
smp<-smp[complete.cases(smp), ]
coordinates(smp)<-~x+y

# fit variogram to sample data
vs <- variogram(sim1~1, data=smp)
m <- fit.variogram(vs, vgm("Exp"))
plot(vs, model = m)

# generate 2 conditional simulations with one core processor
one <- krige(formula = sim1~1, locations = smp, newdata = srgr, model = m,nmax=12,nsim=2)

# plot simulation 1 and 2: statistics (min, max) are ok, simulations are also ok.
spplot(one["sim1"], main = "conditional simulation")
spplot(one["sim2"], main = "conditional simulation")

# generate 2 conditional with parallel processing
library("parallel")
no_cores<-detectCores()
cl<-makeCluster(no_cores)
parts <- split(x = 1:length(srgr), f = 1:no_cores)
clusterExport(cl = cl, varlist = c("smp", "srgr", "parts","m"), envir = .GlobalEnv)
clusterEvalQ(cl = cl, expr = c(library('sp'), library('gstat')))
par <- parLapply(cl = cl, X = 1:no_cores, fun = function(x) krige(formula=sim1~1, locations=smp, model=m, newdata=srgr[parts[[x]],],  nmax=12, nsim=2))
stopCluster(cl)

# merge all parts
mergep <- maptools::spRbind(par[[1]], par[[2]])
mergep <- maptools::spRbind(mergep, par[[3]])
mergep <- maptools::spRbind(mergep, par[[4]])

# create SpatialPixelsDataFrame from mergep
mergep <- SpatialPixelsDataFrame(points = mergep, data = mergep@data)

# plot mergep: statistics (min, max) are ok, but simulated maps show "vertical lines". i don't understand why.
spplot(mergep[1], main = "conditional simulation")
spplot(mergep[2], main = "conditional simulation")


texts: Setting input layers as PostGIS tables when building a model?

texts:
I am trying to build a QGIS 3.8 model where the same inputs are used each time and only the clipping layer extent changes.

I can do this using shapefiles and editing the script to use the same input layer directly but now the job I have been tasked is to do the same thing, only using layers directly from a PostGIS database.

I have seen the following code from another question here which allows me to load PostGIS layers through the Python console but cannot get the layer to act as the input layer itself by changing the Default value in the processing script to the variable uri.

tablename2 = "Sewers_2015 - Outfall"
geometrycol = "geom"

from qgis.core import QgsVectorLayer, QgsDataSourceUri
uri = QgsDataSourceUri()
uri.setConnection("10.0.0.6", "5432", "Outfall", "David", "Pass")
uri.setDataSource ("Outfall", tablename, geometrycol)
vlayer=QgsVectorLayer (uri .uri(False), tablename, "postgres")
QgsProject.instance().addMapLayer(vlayer)

When adding this code to the top of my Clip script it only loads in the PostGIS layer to my workspace but doesnt come up as the default inpput to be clipped.

Here is my full script

user2= 'david.laverty'
tablename = "Sewers_2015 - Outfall"
geometrycol = "geom"
import os
from os import path
from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterVectorLayer
from qgis.core import QgsProcessingParameterFeatureSink
from qgis.core import QgsProject
import processing

from qgis.core import QgsVectorLayer, QgsDataSourceUri
uri = QgsDataSourceUri()
uri.setConnection("localhost", "5432", "Data", "David", "Pass")
uri.setDataSource ("Schema", tablename, geometrycol)
vlayer=QgsVectorLayer (uri .uri(False), tablename, "postgres")
QgsProject.instance().addMapLayer(vlayer)

class SaveToFileTest1(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterVectorLayer('siteboundary', 'site boundary', types=[QgsProcessing.TypeVectorPolygon], defaultValue=None))
        self.addParameter(QgsProcessingParameterVectorLayer('vectors2016', 'Outfall_2016', types=[QgsProcessing.TypeVectorAnyGeometry], defaultValue=uri))
        self.addParameter(QgsProcessingParameterFeatureSink('Clipped', 'clipped', type=QgsProcessing.TypeVectorAnyGeometry, createByDefault=True, defaultValue=os.path.join(r'C:\Users', user2, r'Desktop\pyGIS Testing\clip_test.shp')))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(1, model_feedback)
        results = {}
        outputs = {}

        # Clip
        alg_params = {
            'INPUT': parameters['outfall2016'],
            'OVERLAY': parameters['siteboundary'],
            'OUTPUT': parameters['Clipped']
        }
        outputs['Clip'] = processing.run('native:clip', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Clipped'] = outputs['Clip']['OUTPUT']
        return results

    def name(self):
        return 'save to file test1'

    def displayName(self):
        return 'save to file test1'

    def group(self):
        return 'DL_testing'

    def groupId(self):
        return 'DL_testing'

    def createInstance(self):
        return SaveToFileTest1()

So my question is can a PostGIS layer be set as the Inpupt for a QGIS model using a Python script or am I attempting something that cannot be done as of yet in QGIS.

I'm running QGIS 3.8 on Windows 10 with a PostgreSQL 9.6 install running locally on my machine.


texts: Adding dashed line in QGIS print composer?

texts:
I want to display two maps in the print composer: a main one, and another showing a zoomed in section of the main map. I've drawn a rectangle over the original region and want to add two dashed lines connecting the rectangle to the zoomed in section of the map, but I can't figure out how.

I tried using an arrow without an arrowhead, but there isn't an option to make it dashed. I've also tried drawing a flat dashed rectangle, which works but is quite fiddly.

Does anyone know of an easier way to perform this simple task?


texts: SDE connections to geodatabase as schema owners not working properly

texts:
An enterprise GDB that was working fine seems to have gone wrong somewhere. When I connect to it as the system admin from ArcGIS Pro, I can see all the spatial data that is stored in the SQL Server database properly (vector layer as vector, raster layer as raster). However, when I connect to it as a schema owner (no matter which one), data cannot be seen properly. A bit convoluted to explain, but basically I would say that, when connected as a schema owner:

Raster tables are seen as vector layers (a polygon with the extent of the raster).
Vector layers seem to be recognized properly though (and displayed).
I can't load raster files to the GDB.
If I right click on the connection as a schema owner I can see that "Geodatabase Connection Properties" is greyed out (when I connect as SA, it is NOT greyed out).

I highlight this thing about the "Geodatabase Connection Properties" because I have another Enterprise GDB that doesn't have any problems, and when doing right click as a schema owner, "Geodatabase Connection Properties" is NOT greyed out.
User login name is the same as schema.
-> SQL Server 11.0.2100.60
-> SQL Server 2014 Management Studio
-> ArcGIS Pro 2.6.1


texts: Changing data source of table using ArcObjects?

texts:
Cross post from: https://geonet.esri.com/thread/125010

Does anyone know how to change the data source of a table using arcobjects? I know it can be done with arcpy, but I'm using C#. I've managed to change the source of all the layers in the mxd, but I need to do tables too.

I've come up with a solution by adding and removing from the table collection.

Is this the best way to do it?

        IMap map = ...;
        IWorkspace = ...;

        ITableCollection tableCollection = (ITableCollection)map;
        int tableCount = tableCollection.TableCount;

        // iterate over a copy of the tables so we can add and remove
        List<ITable> tables = new List<ITable>();
        for (int i = 0; i < tableCount; i++)
        {
            tables.Add(tableCollection.get_Table(i));
        }

        foreach (ITable oldTable in tables)
        {
            IDataset dataset = (IDataset)oldTable;
            IDatasetName datasetName = (IDatasetName)dataset.FullName;
            string workspaceName = datasetName.WorkspaceName.PathName;

            // get the new table
            ITable newTable = ((IFeatureWorkspace)workspace).OpenTable(datasetName.Name);

            // replace the table in the TOC
            tableCollection.RemoveTable(oldTable);
            tableCollection.AddTable(newTable);

            // fire change event to change table
            IMapAdmin2 mapAdmin2 = (IMapAdmin2)map;
            mapAdmin2.FireChangeTable(oldTable, newTable);
        }

Note: this is simplified - I did check the table was in the new workspace before opening it.


texts: Automatically generating raster maps (from OSM) for areas bounded by polygons?

texts:
In order to create walking lists for a volunteer campaign I have used qgis to draw polygons over groups of addresses.  These are irregular areas, so I've used postgis's ST_ENVELOPE function to generate rectangular bounding boxes for these areas.

I am looking for way to automatically generate fixed-size raster maps for each of these rectangular bounding boxes.  That is, as output I want to have files named something like map-for-bb-1.png, map-for-bb-2.png, etc.

Doing this manually seems to be quite time consuming (zoom to an appropriate level of detail on the OSM map, save as a raster image, export the rectangular bounding box as a shapefile, use that as a clipping mask for the raster layer, etc...)

How would I go about automating this?  It seems like some of the steps could be performed on the command line (e.g., using pgsql2shp to export each bounding box to an individual shapefile), but I'm not sure the best way to accomplish some of the remaining steps (e.g., downloading OSM raster data for the corresponding rectangle).


texts: SVG symbols are not available on other computers

texts: I created some layers with SVG symbols. Those symbols are saved on the server, as it is not possible to store it in the database (postgis). Every computer can access to the server. But if I open the project on other computers - instead of the SVG symbols there are question tags.


texts: Creating DEM from points, contours and drainage network using Kriging interpolation In ArcGIS Spatial Analyst?

texts:
How to create DEM from points, contours, and drainage network using a kriging interpolation?

The scanned topographic map, with 300 dpi resolution,was georeferenced and manually vectorized using ESRI software package (ArcGIS) and specific extensions. GIS layers containing contour lines, elevation points, and the hydrographic stream network were stored in ESRI shapefile format.

From altimetry data, the digital elevation model (DEM) in GRID format was produced by an advanced geostatistical interpolation techniques such as kriging (Burrough and McDonnell 1998) with a resolution of 10 m.


texts: How to generate random GeoJSON polygons (latitude, longitude)

texts:
I am looking how to generate polygons in GeoJSON format with random locations (latitude, longitude), if I generate totally random set of points it will not produce valid polygons.

So my question is there a tool or algorithm that allows to generate valid polygons of specific size (radius?) inside of bounding box.

UPDATE#1

Here is more details to clarify my question:
Desired polygons should consist of 30-60 points each and represent something similar to square.
All generated polygons should not intersect with each other (no inner polygons either).
All generated polygons should be placed within specified bounding box (border)


texts: How to calculate distance between two locations using stored procedure in MySQL?

texts:
I have a database with locations with their lat and long coordinates. The area of locations is quite small with lots of them quite near each other.

I'm using this stored procedure in my project. http://www.artfulsoftware.com/infotree/queries.php?&bw=1680#109

I'm having troubles with accuracy. I need it to return more accurate results (in metres) instead of just 0 distance.

Can it be modified so I get better results?


texts: Strange slow behavior of a very simple query utilizing single PointField field

texts:
Infrastructure on AWS RDS:

PostgreSQL= 9.6.1

POSTGIS="2.3.0 r15146" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 2.1.1, released 2
016/07/07" LIBXML="2.9.1" LIBJSON="0.12" RASTER

Problem

this is the following snippet of my django class (it has other fields I skipped)

class User(AbstractUser, TimeStampedModel):

    location = gis_models.PointField(default=Point(x=0, y=0, srid=4326),
                                     blank=True,
                                     srid=4326)
    location_updated = models.DateTimeField(null=True,
                                            blank=True,
                                            editable=False)

When I run a pure select query like

SELECT * FROM Users_user where Users_uder.id=1

I get 77 ms execution on the query. If I remove location field (i.e. all but location field) from the query - I get 3 ms response. If I put it back - again the query time grows ~20 times.

Any ideas?

performance debugged through django-debug-toolbar

DATABASES = {
  'default': {
    'ENGINE': 'django.contrib.gis.db.backends.postgis',
    'NAME': 'dbname',
    'USER': 'dbuser',
    'PASSWORD': 'dbpass',
    'HOST': 'postgres',
    'PORT': '5432',
  }
}


texts: ArcGIS saving blank file while using ArcPy to join shapefile and Excel table?

texts:
I'm working on a project that utilizes ArcPy to join a shapefile and a excel file. This join takes information from the excel file and creates a map that highlights certain areas different shades of a color to display the data in a more visually concise format.

I start with a blank map file, and place the shapefile and the excel file into the file as layers. The join works perfectly with no errors; however, when I go to save the map file as a copy, and then subsequently export a pdf of the output, I get blank files both times.

I call arcpy.RefreshTOC() and arcpy.RefreshActiveView() to attempt to make my changes refresh onto the map, but nothing happens.

Here is my code:

import arcpy

# setup workspace
env.workspace = r"C:\Users\Class\Desktop\MapCreation\Resources"
env.qualifiedFieldNames = False

# start work on blank map
mapFile = r'C:\Users\Class\Desktop\MapCreation\Resources\blankMap.mxd'
mxd = arcpy.mapping.MapDocument(mapFile)

# create a table view from the xls
table = arcpy.mapping.TableView(r'C:\Users\Class\Desktop\MapCreation\Resources\2012tpreg.xls\teens$')

# create a new layer using shapefile
shape_file_layer = arcpy.MakeFeatureLayer_management(r'C:\Users\Class\Desktop\MapCreation\Resources\MSshapefiles\MScounties.shp')

# AddJoin
arcpy.AddJoin_management(shape_file_layer, "CNTY_FIPS", table, "FiPS")

# refresh the TOC and ActiveView in order to propagate the changes
arcpy.RefreshTOC()
arcpy.RefreshActiveView()

# export PDF
mxd.saveACopy(r'C:\Users\Class\Desktop\MapCreation\Resources\newMap.mxd')
arcpy.mapping.ExportToPDF(mxd, "test.pdf")

I should note that I know the join is working correctly, because when I run that piece code inside the ArcGIS environment (rather than in IDLE or Command Prompt) it creates the join as a layer in the map environment.

I have tried a lot of different ideas to try to get this problem fixed and I couldn't figure out anything that helped. If anyone has any ideas on what could be wrong I would love to hear them.


texts: Calculating number of polygons needing to be crossed?

texts:
I have a bunch of polygons (parcels, green areas) and a bunch of wells (point, black triangles).

How do I calculate the number of parcels I need to cross to get to the closest well to each parcel?

Straight line distance from the centroid of the parcel is fine.

I am using ArcGIS Desktop 10.6.


texts:
What I mean is, if I have a map which uses the OSGB36 (or EPSG:27700) projection, and I draw two areas (let's say polygons) on my map which cover the same area on the map, does this mean that the real life polygons they correspond to will have the same area in hectares, regardless of the location in Great Britain?
I'm wanting to know this because I'm working with GeoTIFF files using OSGB36, and I want to know if I can calculate the area of a polygon simply by counting pixels within the polygon and multiplying by the real-life area of one pixel.
As an aside, if OSGB36 does preserve area, presumably WGS84 does not - in which case am I right in thinking that it's easier to use OSGB36 if I need area calculations to be accurate?


texts: Creating a "Data entry table" parameter - Arcpy tool scripts

texts:
I want to create an Arcpy tool script that takes values from a 2-column table placed in the tool itself, which is filled by the user, and stores these values in a real table.
In order to do so, I have created a "Record Set" parameter in the script properties. This RecordSet has a data schema taken from an existing table and is composed by 2 columns (param1, param2) (see image below).

Once I execute the tool, the following form appears:

As you can see, the tool has a table where the user can write any values, entering them directly to the table's rows (red dots in the image).
I have come to the solution of getting these values and storing them to a real table using arcpy:
tb = arcpy.CopyRows_management(in_recordset, "table_path")
And there is no problem with this, because it is done quite well.
However, this solution turns out to not being user friendly: In the first place, the upper box with the string "Script::t" on it can't be modified because when doing so, tool crashes. Secondly, one must click first to the "+" button to create a new row, and then fill the spaces. Additionally, while doing so, some cells turn in white and I have to "battle" a little bit with some clicks in order to have them all filled.
So, my question(s):
What is the upper box with the string "Script::t" intended for? And can it be taken out somehow?
And most important: is this the most apt solution to create the tool I need? Are there other options to insert a table to the tool by using another kind of parameter than "Record Set" parameter?
Hope I have explaied myself.


texts:
I'm not sure I'm doing this right, I've two locations, one a point and one a point within an area. I'm running them through st_clusterwithin, and they both create 2 individual cluster groups (which is fine). Problems start when I try to rejoin the cluster using st_contains. Anybody know what I'm doing wrong, or have I stumbled across a weird bug?

-- postgresql-9.5-postgis-2.2
CREATE TABLE test (name character varying, p geometry(point, 900913), a geometry(Geometry, 900913));
INSERT INTO test values ('a', '010100002031BF0D00713D0AD7F54711C1A4703D1A3D965941', null);
INSERT INTO test values ('b', '010100002031BF0D0002FB34FAC62E11C1FAF7FB8141965941', '010300002031BF0D00010000000D0000009A9999998B3211C114AE47F13D965941295C8FC2863211C1A4703D5A569659418FC2F528D33011C1EC51B84E5A9659417B14AE47B62F11C1B81E851B5A965941295C8FC2BC2F11C18FC2F5F85D96594114AE47E10D2C11C1295C8FE26F965941EC51B81EF32B11C1666666D6699659417B14AE47D12C11C1F6285CFF139659418FC2F528BA2F11C1F6285C7F14965941B81E85EBBE2F11C11F85EB01279659415C8FC2F5DB2F11C1D7A3702D2C965941F6285C8F2F3011C15C8FC235309659419A9999998B3211C114AE47F13D965941');

select name, st_contains(ST_CollectionHomogenize(c.g), u.p)
    from test u,
        (select unnest(st_clusterwithin(coalesce(u2.a, st_expand(u2.p, 350)), 250)) g
        from test u2) c order by name;

  name | st_contains
------+-------------
 a    | t
 a    | f
 b    | t
 b    | t
(4 rows)

Notice that the temporary table below is identical to the inner select above:

create temporary table test_g as
  select unnest(st_clusterwithin(coalesce(a, st_expand(p, 350)), 250)) g
  from test;
select name, st_contains(ST_CollectionHomogenize(c.g), u.p)
    from test u, test_g c order by name;
osm_gb-#  name | st_contains
------+-------------
 a    | f
 a    | f
 b    | t
 b    | t
(4 rows)

Neither of the outputs appear to be correct, I'd expect one 'a' row to be true for it's cluster and false for the other, and 'b' likewise.


texts: Removing visible lines between merged Sentinel-2 L2A tiles

texts:
My project area is covered by 4 different S2A MSIL2A tile but when I merge those different tiles on QGIS I can see visible lines between them as shown below.

TILE 1: S2A_MSIL2A_20201115T080211_N0214_R035_T38SKJ_20201115T120134
TILE 1: S2A_MSIL2A_20201115T080211_N0214_R035_T38SKH_20201115T120134
TILE 3: S2A_MSIL2A_20201115T080211_N0214_R035_T38SLJ_20201115T120134
TILE 4: S2A_MSIL2A_20201115T080211_N0214_R035_T38SLH_20201115T120134

How can I remove those lines? With my limited knowledge I couldn't find a way to do it, do I need to apply any correction, if so there is any guidelines I can find on Internet?


texts: Opening a featureclass not being the owner of the data in SDE

texts:
How can you open a featureclass if you are not the schema owner but have been given privs to select or view the data?

So given a name of a featureclass, say "STREETS" - how can I open this if I am using an SDE connection with username "VIEWER" (given privs to see "EDITOR" owned data) but the owner of "STREETS" is "EDITOR"? - using the following code:

IFeatureWorkspace.OpenFeatureClass(string name)

This will work if you are the owner, you just put:

IFeatureWorkspace.OpenFeatureClass("STREETS")

But if you are not the owner then you have to fully qualify it, a la:

IFeatureWorkspace.OpenFeatureClass("EDITOR.STREETS")

Is there anyway round this? Bearing in mind my users will not have access to an sde connection for the data owner.


texts: How to list..only the Field names for a group of files in a folder?

texts: I am having...for eg, 27 files in a folder (may be a shp or tab). Is it possible to list the Field_Names (ID, Type, Date, Allotment, etc..) alone  of all those 27 files in a csv format or xls format?


texts:
How can I split all the lines in a feature class where they intersect with a second line feature.

I have a roads feature and a rail feature. I want the lines in the rail feature to split where ever it crosses a road. I'm looking to do it all at once with an overlay analysis tool. Like the Split Line at Point tool, but using lines as the input. Too many to do it manually.


texts: ArcGIS Repair geometry doesn't solve "Incomplete void poly" error

texts:
I'm working with 8 shapefiles at approximately 63 mb in size total and I am trying to use ArcGIS/arcpy's Union to combine them into one single dataset.
Each time I run this, I get the error:

ERROR 999999: Something unexpected caused the tool to fail. Please refer to the error help for potential solutions, or contact Esri Technical Support http://esriurl.com/support.
Invalid Topology [Incomplete void poly.]
Failed to execute (Union).

On this page, as well as on ESRI's website, it is suggested to use arcgis Repair Geometry to fix this issue. I have tried to run check geometry and repair geometry both with ESRI's validation method and with OGC's validation method. OGC found some errors that it reported solved.
This did however not solve the error.
I have also tried to use v.clean in QGIS 3.2, but this messed up some of the shapefiles geometry (I have a backup, luckily).
Is there some other method i could try?


texts: How to create a filter expression from selection QGIS 3.6.3

texts:
I need to create a filter expression to filter only the selected polygons in a layer. Yes I know I can save the selected elements to a new layer but this solution is not suitable for me.

I managed to create a virtual field with is_selected() that return true only on selected elements, but I cannot use this virtual field in filter composer window.

I really need a "filter from selection" workaround!


texts: Clip raster by mask layer returns error "wrong number of parameters"

texts:
Using Python for QGIS 2.14 Essen i try to clip a raster by a Polygon, as indicated here:
http://docs.qgis.org/2.6/de/docs/user_manual/processing_algs/gdalogr/gdal_extraction/cliprasterbymasklayer.html

with this code:

processing.runalg('gdalogr:cliprasterbymasklayer', input, mask, no_data,     alpha_band, keep_resolution, extra, output)

My code read as following:

processing.runalg('gdalogr:cliprasterbymasklayer', inputlayer, maskshape, "none", False, False, "", outputraster)

Exactly as derived from tutorials. However, this is the error message returned and I can't figure why:
Error: Wrong number of parameters

I´ve adapted the code to the 16 parameters of the new version, although I´m not sure if correct

processing.runalg('gdalogr:cliprasterbymasklayer', inputlayer, maskshape, "", False, False, False, 0, 0, 0, 0, False, 0, False, "",  "outputraster")

The error message says still the same.


texts: pyGRASS on MacOS: unable to import pygrass.raster

texts:
I am running GRASS 7.3 on Mac OS 10.11.5 El Capitan. I'm trying to figure out how to get GRASS GIS running via python scripting (pygrass). I'd like to be able to access raster data programatically. For this I want to use either RasterRow or RasterSegment classes from pygrass.

When I try from grass.pygrass import raster it fails with the following back trace:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/user/Desktop/Mac_pyGRASS.py", line 52, in <module>
    from grass.pygrass import raster
  File "/Applications/GRASS-7.3.app/Contents/MacOS/etc/python/grass/pygrass/raster/__init__.py", line 14, in <module>
    import grass.lib.gis as libgis
  File "/Applications/GRASS-7.3.app/Contents/MacOS/etc/python/grass/lib/gis.py", line 23, in <module>
    _libs["grass_gis.7.3.svn"] = load_library("grass_gis.7.3.svn")
  File "/Applications/GRASS-7.3.app/Contents/MacOS/etc/python/grass/lib/ctypes_loader.py", line 64, in load_library
    raise ImportError("%s not found." % libname)
ImportError: grass_gis.7.3.svn not found.

I am able to import other pygrass modules without error.

The relevant portion of my script is as follows:

#!/usr/bin/env python
import os, sys
import subprocess

# path to the GRASS GIS launch script (MacOS)
grass7bin = '/Applications/GRASS-7.3.app/Contents/MacOS/grass73'

## SPECIFY LOCATION AND MAPSET
## both location and mapset must already exist
location = 'mylocation'
mapset = 'PERMANENT'

# query GRASS 7 for location of GISBASE
startcmd = [grass7bin, '--config', 'path']

p = subprocess.Popen(startcmd, shell=False,
                    stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, err = p.communicate()
if p.returncode != 0:
   print >>sys.stderr, "ERROR: Cannot find GRASS GIS 7 start script (%s)" % startcmd
   sys.exit(-1)
gisbase = out.strip('\n\r')

## SPECIFY GRASS ENVIRONMENT VARIABLES
# Set GISBASE environment variable
os.environ['GISBASE'] = gisbase ##  = '/Applications/GRASS-7.3.app/Contents/MacOS'
# the following not needed with trunk
os.environ['PATH'] += os.pathsep + os.path.join(gisbase, 'extrabin')
# add path to GRASS addons
home = os.path.expanduser("~")
os.environ['PATH'] += os.pathsep + os.path.join(home, '.grass7', 'addons', 'scripts')

# define GRASS-Python environment
gpydir = os.path.join(gisbase, "etc", "python")
sys.path.append(gpydir)

# Define GRASS DATABASE
# add your path to grassdata (GRASS GIS database) directory
gisdb = os.path.join(os.path.expanduser("~"), "grassdata")
# Set GISDBASE environment variable
os.environ['GISDBASE'] = gisdb

# Import GRASS Python bindings (see also pygrass)
import grass
import grass.script as gscript
import grass.script.setup as gsetup
import grass.pygrass
from grass.pygrass import raster

from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r

## launch session
gsetup.init(gisbase,
           gisdb, location, mapset)

# say hello
grass.message('--- GRASS GIS 7: Current GRASS GIS 7 environment:')
print grass.gisenv()

# exit
sys.exit(0)


texts: Integrating Vector/Raster in Google Cloud Storage in QGIS

texts:
I was looking at the new feature of QGIS 3.2 which can integrate between Google Cloud Storage and QGIS. I'd like to save all >20 GB GIS Dataset into Cloud rather than in my local disk. Can I do that?
What should I fill on 'Bucket/container' and 'Object Key' sections?


texts: Moving map-embedded annotation to a different annotation group

texts:
I'm working with an annotation group that was generated by converting labels to annotation and storing it in the map instead of a database.  After this conversion, the scale of the map was changed, and a number of edits were made to the annotations.

I needed to duplicate a few of the annotations to indicate parts of some multi-part polygons that originally only had one label for each part.  When I copied and pasted an existing label in the annotation group, however, I made the mistake of placing it in the Default annotation group instead of the one I was editing.

Since the Default group has a reference scale fixed to the data frame, the pasted annotations are a much different size than the original annotations.  Is there a way to move annotation from one annotation group to another?


texts: (How) can I write string data to a raster with Python/GDAL?

texts:
In Python, I have a NumPy array of string values (categories of some variable) that I'd like to write to a raster file (.tif or .img) using GDAL.

The impression I get is that it's not possible to write a string value into the grid cell of a raster file, meaning I'd need to create an integer lookup table and write the integers into the grid cells instead. Then I'd need to write the character data into the column of the raster's attribute table so that you could know what string values the integers represent.

Is my reasoning correct? If so, I'd like to know the basic process for building a raster attribute table using Python, GDAL. I haven't had much success in the past, and, as I understand it, not all GDAL raster formats are able to have attribute tables at all anyway.


texts: Cannot locate GRASS plugin in QGIS 2.14.1 on Linux

texts:
I am trying to load the GRASS plugin in QGIS but it is not shown in the plugin manager (no matter if I turn on experimental and deprecated plugins). Under processing / options / providers GRASS GIS 7 commands and GRASS commands are visible and activated.

I am running Linux Mint 17.3 Rosa. QGIS 2.14.1 (after upgrading from 2.8. where I had the same problem), GRASS 7, and the qgis-plugin-grass (version 1:2.14.1+20trusty) are installed.

I have read this post How to activate GRASS GIS 7 plugin in QGIS 2.8?. There it is said that QGIS 2.8 is not compatible with GRASS 7. Therefore I installed QGIS 2.14, which should be compatible (http://www.qgis.org/en/site/forusers/visualchangelog212/index.html?highlight=grass#feature-update-of-the-grass-plugin).

Some more information.
I get this errors in the QGIS plugin log:

Failed to load /usr/lib/qgis/plugins/libgrassplugin6.so (Reason: Cannot load library /usr/lib/qgis/plugins/libgrassplugin6.so: (libgrass_gis.so: cannot open shared object file: No such file or directory))

Failed to load /usr/lib/qgis/plugins/libgrassprovider6.so (Reason: Cannot load library /usr/lib/qgis/plugins/libgrassprovider6.so: (libgrass_gis.so: cannot open shared object file: No such file or directory))

Failed to load /usr/lib/qgis/plugins/libgrassrasterprovider6.so (Reason: Cannot load library /usr/lib/qgis/plugins/libgrassrasterprovider6.so: (libgrass_gis.so: cannot open shared object file: No such file or directory))

So, it seems QGIS looks for Grass 6 plugins while Grass 7 is installed. I can set LD_LIBRARY_PATH, but then the Plugin (6) does not match Grass (7). In /usr/lib/qgis/plugins only the lib...6.so are available. How can I get the 7?


texts: USGS topo symbols for light-duty and unimproved roads

texts:
Mapping in QGIS 3.14, I want to use the same symbol as the local historical USGS 24K quads use for unimproved road and light-duty road.
An example can be seen here (the black continuous road and the black dashed roads):

For an unimproved road, that is an interrupted double line (like a string of equals signs) with no fill: the map background shows through. I haven't found anything close in QGIS.
For a light-duty road, that is a continuous double line (one long, twisty equals sign) with no fill. QGIS has something similar, except that it has white fill, brightly painting over the map behind it.

Am I missing something that is already in QGIS 3.14?

Can I obtain these two effects from pieces that already are in QGIS 3.14?

Is there an appropriate symbol library that I can obtain for free or money?


texts: Reading ASCII Grid file cell value for given long/lat coordinates

texts:
I am trying to write a Python script which takes a set of long/lat coordinates as input and returns the  cell value corresponding to these coordinates

The ASCII Grid file is from PVGIS and contains solar radiation data.

Metadata:

Format: ESRI ascii grid
Map projection: geographic (latitude/longitude), ellipsoid WGS84

The header of the .asc file has this format:

ncols 5201
nrows 4001
xllcorner -65
yllcorner -35
cellsize 0.025
NODATA_value -9999

Is there a nice way of doing this without ArcGIS ?

I am looking into GeoPandas but have not yet found a way to solve this.

I have done some testing, initially with gdal but gave up on that for several reasons, so decided to make my own script inspired by world2pixel(). See my code below.

So now I have a couple of new questions/issues:

1.
If I input in my function readCellvalue() for example latitude 13.7 and longitude 60.7 my cellY value will become negative, which does not work very well. I am confused.

Should I simply put an abs() on it or is the problem somewhere else?

Can someone help me track this issue?

2.
Could someone clarify "ulY = 0.0 + 6 * 50.0" as suggested by @user1269942

My code:

`

import numpy as np

data = open('filename.asc')

def readCellvalue(lat, long):
   xllCorner = -65.0
   yllCorner = -35.0
   cellsize = 0.025

   lines = data.readlines()

   cellX = int((lat - xllcorner) / cellsize)
   cellY = int((yllCorner - long) / cellsize)

   line = lines[cellY]

   # line is being evaluated/separated at every " "
   myArray = np.fromstring(line, dtype = float, sep = ' ')

   # picks out the right "column" from my array
   cellValue = myArray[cellX]

   if cellValue == -9999.0:
       print("This cell is at sea ! ! ! ")
   elif cellValue != -9999.0:
       print("Cell value found: " + str(cellValue))

`


texts: Batch Dissolving FC into New FCs based on Field Name using ArcPy?

texts:
I am trying to take a FC, export individual FCs based on fields in the original FC, and have the new FCs dissolved.

I found a Python script that works and modified it for my use - see below.

# import system modules
import arcpy
from arcpy import env
import os

# Set environment settings
env.workspace = "H:/GISData/501/Jason/CFRPM7.gdb"

# Set local variables
inFeatures = "YR_2015/Network_2015"
fields = arcpy.ListFields

x = fields
# Execute Dissolve
arcpy.Dissolve_management
(inFeatures, "H:/GISData/501/Jason/Network/node_dissolve.gdb/fc_ + str(x)", ["ROAD_NAME", "NUM_LANES", "POST_SPEED", "AREA_TYPE", "FAC_TYPE", "UA_NAME", "UA_TYPE", "FUNCLASS", "SIS", "NHS", "PAVED", "TRFC_CALM", "SCHL_ZONE", "COUNTY", "MPO", "FDOTDIST", "JURIS", "US_REP", "FL_REP", "FL_SEN", "CNTY_DIST"],"", "SINGLE_PART", "DISSOLVE_LINES")

The problem I'm running into is that I want the exported FCs to be named based on the field they represent. For example, I want the FC exported from the field "ROAD_NAME" from the YR_2015/Network_2015 FC to be named "ROAD_NAME" (.../node_dissolve.gdb/ROAD_NAME). This does not appear to be working, and instead exports the first listed FC, and as (.../node_dissolve.gdb/fc_).

I am new to Python.


texts: Spatial Query for centroid of one feature falls within another

texts:
I am looking at doing a spatial query where I can find when the center of one feature falls within another. Both geometries are on different layers. For example, when the center (centroid) of an area (source layer) falls within the polygon of another area (target layer). If the centroid of the source layer is within the target layer, the expression should then give out the specific attribute from the target layer.
Example: The orange layer is created by myself. Now I want to compute the centroid of that layer. The centroid would fall within the blue layer, here in the Polygon "BSAltona ...". Now in the blue layer polygon there is an attribute I want to extract (a link to a website).

Right now I am using Reffunction (geomwithin and geomverlaps) however there are to many situations, in which the the orange layer would be intersecting multiple of the blue layers (see Picture 2). In this case none of the reffunctions gives out the fitting attribute.
Right now I am using the following:
[%if(geomwithin('layer', 'attribute') is not NULL, geomwithin('layer', 'attribute'), geomoverlaps('layer', 'attribute'))%]

In this situation it would be most logical in my case to give out only the information of the  blue layer that intersects with the centroid of the orange layer. However is think this is not possible with reffunctions.

Also, I am looking to this as an expression only, as I want to insert this into the atlas function. This is why creating a separate layer using the normal create centroid function would hinder my workflow too much.


texts: Automatically create a polygon around other polygons depending on attributes

texts:
I have 1000+ polygons, each has an attribute "Project_Area".

There are maybe 8 different project areas. Is it possible using built-in Tools or Python to create a new shapefile that encompasses each of the 8 project areas?

Its not the hardest task to do manually, but I shall need to do it many times in the next few weeks.


texts:
I am trying to get a search cursor to go through each row and Multi buffer a feature.

I keep getting a syntax error with no error codeother than Invalid Syntax (traceback to this line of code) so i am assuming its something I am doing in the calling of the multiring buffer.

Error Message:

C:\Python27\ArcGIS10.3\python.exe N:\GIS\Projects\AA_Leith_Hawkins_TestBed\Search_Cursor\Script\SearchCUrsor.py
File "N:\GIS\Projects\AA_Leith_Hawkins_TestBed\Search_Cursor\Script\SearchCUrsor.py", line 19
arcpy.MultipleRingBuffer_analysis('Holdings_Layer', 'ofc', distances, unit, "", "ALL")
    ^
SyntaxError: invalid syntax

Code below

import arcpy
from arcpy import env
arcpy.env.overwriteOutput = True

file_workspace = "r'N:\GIS\Projects\AA_Leith_Hawkins_TestBed\Search_Cursor\Search_Cursor.gdb"
env.workspace = file_workspace
Var_Polygonname = row[0]
Holdings = r'N:\GIS\Projects\AA_Leith_Hawkins_TestBed\Search_Cursor\Search_Cursor.gdb\Data\Holdings'
ofc = r'N:\GIS\Projects\AA_Leith_Hawkins_TestBed\Search_Cursor\Search_Cursor.gdb'
distances = [1000, 4000]
unit = "Meters"
arcpy.MakeFeatureLayer_management(Holdings,'Holdings_Layer')

with arcpy.da.SearchCursor(Holdings, ['Holding_Reference_Number'])as Holdings_Ref_cursor:
for x in Holdings_Ref_cursor:
    print x[0] + Var_Polygonname

    arcpy.SelectLayerByAttribute_management('Holdings_Layer', 'NEW_SELECTION', '"Holding_Reference_Number" = \'Var_Polygonname\''
    arcpy.MultipleRingBuffer_analysis('Holdings_Layer', 'ofc', distances, unit, "", "ALL")
   #arcpy.Buffer_analysis("Holdings_Layer", ofc, var_Buffer, "FULL", "ROUND", "ALL", "")
    print'Buffer complete'


texts: Creating series of points along multiple polylines at once in ArcMap?

texts:
I have a series of polylines that represent linear vegetation (shelterbelts, etc). I need to convert these to a series of points equally spaced in order to add them to a different shapefile that represents point vegetation (individual trees). I should also mention that I have about 300 polylines to do this to, so would like to be able to do them all at once. The polylines vary in length, but regardless of length, I want to place points at the start at end of each line, and then every 10 m.

I would imagine that this will require python script. Of which I am not familiar with how to write.

Does anyone know how I can do this in ArcMap 10?

I looked at a similar question previously posted for someone who wanted to create 6 points on 70,000 lines of various length and found this python script:

import arcpy
points = []
  for row in arcpy.da.SearchCursor("YourInputLinesFileHere", ["SHAPE@"]):
    length = row[0].length
    for i in range(0,6):
      point = row[0].positionAlongLine(length/5*i)
      points.append(point)
      arcpy.CopyFeatures_management(points, 'YourOutputPointsFileHere')

So I know it must be possible with python but am not familiar with it enough to manipulate this.


texts: Making bounding box for handling label collision on a map

texts:
In my application, I have to make an arrangement so that overlapping labels for roads/streets are handled correctly.

So I have decided to determine bounding box of each label and then check them for collisions. But I am facing an issue.

What I have with me is the lat-long coordinates of the label aligned along the road.So I have lat-long centre of the bounding box. I also have height and width of the label in meters as well as screen dots. But I am not able to decide what units the rectangle sides should be in to have correct collision detection.

Should rectangle be made with degree measurements or with meters measurements or in terms of screen dots?


texts: Equal area projections with different area results in ArcGIS for Desktop?

texts:
Using ArcGIS, I have to calculate areas of huge polygons all over the Pacific Oceans betwwen 20N and 30S. I have been trying different Equal-areas projections (Cylindrical and sinusoidal for example) and found different results for the area claculation.

For example, one polygon is having a surface of 1'604'176 km2 with Cylindrical and a surface of 1'605'483 km2 with sinusoidal...

So my questions are:

Am i correct to use Equal area projections for calculating areas? Should i be using a different type of projections (equidistance...)?

How can you have different results between the projections and which one should be consider best?

At last, using an equal area projection, I discover that if I break down a big polygon into 6 smaller polygons and i compare the sum of the areas of the 6 smaller polygons with the area of the big polygon, i have a difference between the values (difference of 60Km2 for a polygon of 1'600'000Km2). How is it possible?


texts: How do I create line vectors with a new ID from points in QGIS?

texts:
I have a list of stations from which I would like to create a network of lines (the one created with the Delaunay triangulation would be fine) identified by a particular ID (for example, "AB" for the line from station A to station B) to which associate a list of measurements I already have (distance etc.).

I tried to use polygons to lines on the Delaunay triangulation layer, but the new lines vector has exactly the same table of it and I can't find a way to discriminate all the lines with a new ID...

This is the main table to which I would like to add all the new lines with the ID listed in the first column:

"Fst" and "Matrice di distanza_Distance" are the measurements I need to associate to each line.

These are the tables for the Delaunay triangulation and the Polygons to lines layers, respectively:


texts:
I have a shapefile with UK districts and unitary authorities. In the past, using commuting patterns, I worked out a new geography which is close to travel to work areas.

I now would like to create a new shapefile where UK districts and Unitary authorities are merged to create my new map of travel to work areas.
I know it is possible to merge features manually by clicking the "merge selected features" icon. Only it is taking me ages...

Is there a way of programming this?

Is there a way to write some syntax in order to merge features?

That would help a lot, because I will be able to write the program and while I am running it doing something else. Also I would have the syntax saved somewhere.

I am using QGIS version 1.8.0.


texts: Clipping in PyQGIS2 working but error in PyQGIS3 (argument 1 has unexpected type 'str')

texts:
I am trying to clip streets by a layer with municipalities in them. I was trying to run this script in QGIS3 but obtain the error hereunder.

  Traceback (most recent call last):   File
  "C:\OSGEO4~1\apps\Python37\lib\code.py", line 90, in runcode
      exec(code, self.locals)   File "", line 1, in    File "", line 11, in    File
  "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\tools\general.py",
  line 105, in run
      return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)   File
  "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\core\Processing.py",
  line 130, in runAlgorithm
      context = dataobjects.createContext(feedback)   File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\tools\dataobjects.py",
  line 72, in createContext
      context.setFeedback(feedback) TypeError: QgsProcessingContext.setFeedback(): argument 1 has unexpected type
  'str'

I tried to run the script in QGIS2 and it works perfectly (changing run to runalg off course).

I checked the documentation of clipping in PyQGIS3 but I can't figure out what I am doing wrong. I hope someone has some insights for me.

import processing

streets = QgsVectorLayer(r"C:\Users\Me\Desktop\Test\\streets_fix.shp", "streets", "ogr")
if not streets.isValid():
  print("Layer failed to load!")

municipalities= QgsVectorLayer(r"C:\Users\Me\Desktop\Test\municipalities.shp", "municipalities", "ogr")
if not municipalities.isValid():
  print("Layer failed to load!")

processing.run("qgis:clip", municipalities, streets, r"C:\Users\Me\Desktop\Test\clipttest1.shp")


texts: Is there a way to link ArcGIS to MS Access and use that DB?

texts:
I have an IT background but I am new to GIS. I recently recieved training using ARCMAP.

I have about 5000 farmers that I need to put in a database and access the information using GIS. The vision is to have a map and when I click on a particular plot of land I need to see the related farming information (Owner, productivity etc)

My question is - could I create a database using ARCGIS or is there a way to link ARCGIS to MSAccess and use that DB. I need someone to point me in the right direction. The resulting Map/Database should be Easy to use and update. Reports amay also neeed to be run on a monthly basis.


texts: Shader is not working correctly / set(material.[]) is buggy

texts:
Concider this first Example:
My CGA File:
const dryWallCol            = "Walltexture1/Drywall4.png"

Lot -->                 primitiveCube(20,20,20)

                        setupProjection( 0 , sayit , 37.5 , 25)
                        projectUV(0)
                        set(material.colormap, dryWallCol)

My Texture is mostly white, but in the viewport and later on export to Unreal the texture is a dark Gray. Why?

Concider the second Example:
MyCGAFile:
const leer              = ""
const sayit             = "scope.xy"
const dryWallCol            = "Walltexture1/Drywall4.png"
const dryWallNor            = "Walltexture1/Drywall4Normal.png"
const dryWallHei            = "Walltexture1/Drywall4Height.png"

Lot -->                 primitiveCube(20,20,20)

                        setupProjection( 0 , sayit , 37.5 , 25)
                        projectUV(0)
                        set(material.colormap, dryWallCol)

                        setupProjection( 1 , sayit , 37.5 , 25)
                        projectUV(1)
                        set(material.bumpmap, dryWallHei)

                        setupProjection( 2 , sayit , 37.5 , 25)
                        projectUV(2)
                        set(material.dirtmap, leer)

                        setupProjection( 3 , sayit , 37.5 , 25)
                        projectUV(3)
                        set(material.specularmap, leer)

                        setupProjection( 4 , sayit , 37.5 , 25)
                        projectUV(4)
                        set(material.opacitymap, leer)

                        setupProjection( 5 , sayit , 37.5 , 25)
                        projectUV(5)
                        set(material.normalmap, dryWallNor)

                        setupProjection( 6 , sayit , 37.5 , 25)
                        projectUV(6)
                        set(material.emissivemap, leer)

                        setupProjection( 7 , sayit , 37.5 , 25)
                        projectUV(7)
                        set(material.occlusionmap, leer)

                        setupProjection( 8 , sayit , 37.5 , 25)
                        projectUV(8)
                        set(material.roughnessmap, leer)

                        setupProjection( 9 , sayit , 37.5 , 25)
                        projectUV(9)
                        set(material.metallicmap, leer)

I set every other material to empty, so I can check if they are the problem - they are not. I also checked the Textures several times, the color profiles are allright etc. but as soon as I try to use either of them, the builtin UnknownTexture.png plasters over everything. Why?

Concider my third example:
imagine the code from before, but right after the primitive cube there is
set(material.shader, "CityEnginePBRShader")

This only solves the problem in the first example if you set the bumpmap.
But normalmap still can't be used because the builtin UnknownTexture.png plasters over everything again.


texts: r openproj to plot in open street map in a different projection

texts:
I am very new to geographic data, so I am sorry for the stupidity of the question.

I need to plot some data that are in this projection

EUREF_FIN_TM35FIN
Projection: Transverse_Mercator
False_Easting: 500000,000000
False_Northing: 0,000000
Central_Meridian: 27,000000
Scale_Factor: 0,999600
Latitude_Of_Origin: 0,000000
Linear Unit: Meter
GCS_EUREF_FIN
Datum: D_ETRS_1989

in R in openstreetmap.
I do not understand how the openproj works and looking at the guide here http://cran.r-project.org/web/packages/OpenStreetMap/OpenStreetMap.pdf did not help.

What I came up with is this

a<-min(min(samples$lon))#, min(samples_all$lon))
b<-min(min(samples$lat))#, min(samples_all$lat))
c<-max(max(samples$lon))#, max(samples_all$lon))
d<-max(max(samples$lat))#, max(samples_all$lat))

map = openmap(c(lat= d,   lon= a),
              c(lat= b,   lon= c))
              #              ,minNumTiles=9,type=nm[i])
plot(map)

map_fin <- openproj(map.in, projection =
                      "+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs")
plot(map_fin)

that of course does not work, can you please help and maybe tell me where to find information about projection in R and openproj? I struggled a lot but could not find any.

Also how does this work with ggmap?


texts: Why does GDAL Translate not create any output in batch mode?

texts:
Has anyone tried using the Batch Mode function in GDAL translate - qgis2?  I can process files individually ok, in my case MODISA HDF files. But the batch mode does not seem to work.  Here is my set up parameters

After hitting OK, I get a message instantly saying that the operation is complete, but there is nothing in the Output directory or anywhere else.

Also, oddly, when I process files individually, I've got the 'Load into canvas' option ticked, but I get an error saying that this failed.  Its not such a big deal because I can add it ok, but perhaps this is another bug.  Any suggestions welcome.  Thanks!

Forgot to say that I'm using Windows 7, qgis2.


texts: Merging multiple shapefiles into one shapefile

texts:
I have 5 shapefiles in Postgres. i am trying to merge them into one shapefile but each shapefile has a different amount of columns and data types. I am aware the union and union all will not work if there are different amounts of columns and data types. they all have the same SRID

select * into newtable
from(select * from schools
    union all
    select * from childcare
    union all
    select * from infrastructure
    union all
    select * from nursinghomes
    union all
    select * from hospitals
    )t;

I am wondering if there is work around that can merge or smash the shapefiles together similar to the merge tool in arcgis which will do this


texts: How to loop a Select by Attribute to use in a Select by Location?

texts:
I am new to python scripting and I need help.

I am trying to select one row at a time on a shapefile that contains buffer rings around sites we are working on (some rings overlap), and for each site buffer I need to select by location and extract the castrate parcels that intersect.

At the moment the code runs through ok but it outputs the entire cadastre dataset for each site and does not appear to select each row to run in the select by location function. What am i doing wrong?

arcpy.MakeFeatureLayer_management(Cadastre,'cadastre_lyr')
arcpy.MakeFeatureLayer_management(Buffer_5km, 'buflayer')

field = "SiteName"
cursor = arcpy.SearchCursor('buflayer')
for row in cursor:
    out = (row.getValue(field))
    arcpy.SelectLayerByAttribute_management("buflayer", "NEW_SELECTION")
    arcpy.SelectLayerByLocation_management('cadastre_lyr',"INTERSECT", 'buflayer')
    arcpy.CopyFeatures_management('cadastre_lyr', CadastreOut + "Cadastre_" + str(out) + ".shp")


texts:
Is there a way to source a shapefile to a layer file in QGIS much like how you do this in the properties of an ArcGIS layer? (i.e. Go to layer properties, source tab and then just depress the 'Set Dat Source' button)

I've been all over the layer properties in QGIS and I don't see this anywhere...

EDIT: I've received comments explaining that you cannot do this with shapefiles, but if I were to change the name of a shapefile and then open a pre-existing project file which contains this shapefile as a layer, I will get the 'Handle Bad Layers' dialogue, which allows me to navigate and re-source the layer to any shapefile I want to. I find it odd that you can resource in that scenario, but you cannot do it manually through the layer properties dialogue box.


texts: Formatting numeric fields in attribute table in QGIS

texts:
In ArcGIS you have the following window that allows you to visualize your numeric values in your attribute table in a diffrent manner. (Sorry the Image is in German)

In the example numeric values will be rounded to the 3rd decimal number and the alignment will be left. This does not change the data, it is simply represented differently in the attribute table. in ArcGIS you could also set the settings to show the data in scientific style or in fractions
I would love to do this in QGIS as well. To tell QGIS to represent the data with a max of 2 decimal numbers but actually use all decimal numbers when computing.
Is this somehow possible?
I'm using QGIS 3.4.13


texts: Copying tables results in duplication of some tables, feature classes and relationship classes

texts:
ArcGIS Pro 2.2.3.

I have a file geodatabase called "supportSource.gdb" with 3 tables, 2 feature classes and 3 relationship classes. I just want to copy the contents of this file geodatabase to another file geodatabase using python.  Fairly straightforward, or so I thought.

Tables

EnvGENvecAuditStreetSweeping__ATTACH
EnvtblTreeMaintenanceRecords
EnvtblTreeMaintenanceRecords__ATTACH

FeatureClasses

EnvGENvecAuditStreetSweeping
EnvvecTreeMaintenanceAreas

Relationship Classes

EnvGENvecAuditStreetSweeping__ATTACHREL
EnvtblTreeMaintenanceRecords__ATTACHREL
TreeMaintenance

If I try selecting all the tables and copying from the source geodatabase to a new empty geodatabase (called "dest.gdb" in Catalog, I get the following error:

Next, I tried a simple python script to copy the tables only into the destination geodatabase.

import os
from arcpy import env
destgdb = "E:\\Temp\\dest.gdb"
sourcegdb = "E:\\Temp\\supportSource.gdb"
env.workspace = sourcegdb
tableList = arcpy.ListTables()
for table in tableList:
    arcpy.Copy_management(table, os.path.join(destgdb, table))

Not only does this copy the tables, it also copies the feature classes as well as the relationship classes. This is not a problem for me because I need to copy all contents, not just tables. However, what does occur is a duplication of one of the tables, one of the feature classes, and one of the relationship classes. It also renames one of the relationship classes weirdly ("EnvtblTreeMaintenanceRecords__ATTACH__ATTACHREL" instead of "EnvtblTreeMaintenanceRecords__ATTACHREL")

My guess is that one of the tables, EnvtblTreeMaintenance participates in two relationship classes. It is the Destination for the TreeMaintenance relationship class and the Origin for the EnvtblTreeMaintenanceRecords__ATTACHREL relationship class.

This behaviour is not seen in ArcGIS Desktop 10.6.1. When I copy the three tables manually in Catalog in Desktop, they copy across without problems. All feature classes and relationship classes also copy along with the tables. However, when I use the above python code in Desktop, everything seems to copy correctly, but I get the following error in the python console:

  Runtime error  Traceback (most recent call last):   File "",
  line 2, in    File "c:\program files
  (x86)\arcgis\desktop10.6\arcpy\arcpy\management.py", line 4314, in
  Copy     raise e ExecuteError: ERROR 000725: Output Data Element:
  Dataset E:\Temp\dest.gdb\EnvtblTreeMaintenanceRecords__ATTACH already
  exists.

My question: what is the best way to copy the data in the source geodatabase to the destination database in ArcGIS Pro without duplication of tables, feature classes and relationship classes?


texts: Unloading addresses from database to Nominatim/Pelias geocoder?

texts:
I don't works with geocoding before. But I have geocoding task now. We have postgreslq database with addresses data. And we have two ways for geocoding: nominatim and pelias (are placed on different server than address DB).

There is issue:

how unload addresses data from data base for geocoder well?
We have to choose .csv or something more?
maybe exists way of linking DB with geocoder without unloading?


texts:
I'm hoping someone out there might be able to help or suggest some tests to try and resolve some very odd QGIS behaviour of late...

Our core GIS team of 4 users has been testing v2.10 prior to roll-out for the remainder of the organization.

In the past couple of weeks we have had some very odd shapefile corruption things happening and are yet to isolate or identify a potential cause...

These have included:

Editing one line (using the scissor tool to split into two records) and having the resulting records attributes completely replaced by attributes from another geometry no-where near the edited row (spatially or in the attribute dbf)
Line geometries completely disappear while their attributes get associated with an unrelated geometry
Line geometry get partially deleted and replaced with a vertice in infinity (line visually extends beyond the dataset and no matter how far zoomed out you never see the end)

The above is happening in a polyline layer, that we did discover has some multi-line geometries..? Perhaps this is impacting... but never saw this issue with v2.8 and the same layer.
The layer has undergone MANY save-as's in the past week and many verify geometry checks as we thought invalid geometry may have been the cause...

In another separate layer yesterday, to a separate user (me) i created a layer with approx 30 polygons... adding, deleting some and doing some 'avoid intersection' to remove overlaps. Layer looks fine in qgis, saved and also a save-as to a completely new file but when i load into ArcGIS online (and later is ArcMap as a test) some of the overlaps and deleted objects are still there. Close and reopen qgis,  reopen the file and they are definitely not.

The project was created in v2.8 that i was editing in (but had been saved in 2.10).

Does anyone have any ideas/suggestions.

We do still have 2.8 installed as well as 2.10? We have never encountered issues like this when testing other new qgis versions while still running old ones?


texts: Technique using GDAL to overlay / fit one PNG map to directly overlay on another?

texts:
I have one .PNG basemap with known North edge of Latitude = 52.48, a South edge of Latitude =  21.94, an east edge of Longitude = -61.87 and a west edge of Longitude = -129.37.

And I have another PNG image overlay with a known North edge of Latitude = 50.03, a South edge of Latitude =  22.02, an east edge of Longitude = -67.24 and a west edge of Longitude = -124.99.

Is there a way using GDAL or other simple scripting tools to process the second PNG image to be in an output that directly overlays (same x/y dimensions) the basemap?

I have no trouble generating map tiles that overlay each other beautifully - but I seem to be having no luck figuring out the syntax of the seemingly simpler case.

Help in the form of a simple example is much appreciated!

To be clear: the goal is an output .png file that has the exact same x/y dimensions as the basemap, and which can directly overlay on top of it.


texts:
i want to diplay a kml (on a public url) on this map
http://rumsey.geogarage.com/maps/g2104028.html?lat=44.22158376545796&lon=10.634765625&zoom=8

any hint to start?
looking at that page source i see
var map = new geogarage.Maps(document.getElementById("map"), MY_IMAGE, {method:"overlay"});

but i'm not finding anything about geogarage documentation


texts: Correct Incorrect Projection of Imported Shapefile into Excel PowerMap

texts:
Update 2
Below are the steps I re-projected shapefile. For some reason it disappears.
In QGIS double-click on shapefile

Right-click and select Set Layer CRS

Select WGS 84

Click OK

Right-click shapefile and select "Save As"

Create new subdirectory to save new shapefile

Click OK in Save Vector Layer

I uncheck the original shapefile and the new shapefile is blank

Update 1
Still struggling, found some clues.
This website, Power Map for Excel gives simple example and links to the Excel spreadsheet and shapefile. When I plug it in, it looks correct.
I went back to my original shapefile and it is plotting the boroughs of New York City across the world.
In QGIS the shapefile looks like this

But when I import into Excel Power Maps it looks like the following.
How to correct this? They are all points on a map instead of shapes AND they are all in the wrong place!
Staten Island in the North Pole

Bronx is South of South Africa

Queens is in the Indian Ocean

Brooklyn is North of the Philippines

Manhattan is North of Indonesia

Original Question
I am very novice at shapefiles and cartography in general.
Finally found a shapefile (shapefile is under Borough Boundaries & Community Districts called "Borough Boundaries Clipped to Shorelines") that does not give errors when I import into Excel 2016 Power Map
I am not sure how to use Shapefile that I successfully imported.
I googled and searched for YouTube. YouTube only shows animation of completed maps, not a how-to.
And the book High Impact Data Visualization in Excel with Power View, 3D Maps, Get & Transform and Power BI, 2nd Edition (free download at allitebooks) gives one or two sentences, not even a tutorial on how to import for use!
Here are the steps, perhaps you will understand
I opened Excel 2016 Power Map > Clicked on Custom Regions

Click on Import New Set

Select BoroName

Click on Import

Map view defaults to the Indian subcontinent, which I find odd because the shapefile is 5 boroughs of New York

I attempt to add a field from my data whose name is BOROUGH. Entire dataset can be downloaded from NYS Open Data

Then I add Custom Region Set and nothing happens.
What am I supposed to do next?


texts: How can I convert SRTEXT to PROJ4TEXT with QGIS

texts:
I have a SRTEXT like this: GEOGCS["My Data",DATUM["Not specified (.....

I want to get PROJ4TEXT of this STREXT.

Can I get it from QGIS. Because I want to insert it in Postgresql SPATIAL_REF_SYS table. (or I need PROJ4TEXT data to use this projection)


texts: Extracting locatation information from Google Maps and business directories

texts:
Is there a way to mass download/extract (infrastructure) data from either Google Maps or Business directories based on a search query?

I am often asked to create inventories of assets within a local authority such as places of worship, schools or hospitals. Currently I just do a Google Map search and extract each name and address manually before geo-referencing in Mapinfo which is very time consuming. Is there a tool which will download all this information automatically (name, address, postcode, phone number) from Google Maps or 118118, or is copy/pasting the results the only way.


texts: Not able to apply update cursor on OS MasterMap feature classes created using OS Translator II

texts:
I am using OS Translator tool to convert OSMasterMap GZ files to PostgreSQL database. The feature classes created are query layers and cannot be edited in ArcGIS.
Now I am writing a code to add new data to the existing feature classes. For this purpose, I am deleting the old rows in the given extent of new data and inserting the new data.
However, I am getting an error "Create cursor failed" while deleting and I haven't even tried Inserting.

What I have done as an alternative workaround is that :

I created shapefiles of the OS Master Map feature classes locally.
Then I created new feature classes in PostgreSQL with similar fields which is for some reason editable.
I am then trying to insert the data from shapefile to this new editable feature class.

After this I am able to update and insert data on the new feature class when new data comes.

Please can you tell me if I can directly update and insert the feature classes created by OS Translator so that I can skip these steps in the middle.

Here is my code:

import arcpy
import datetime

temp_workspace = arcpy.GetParameterAsText(0)
sde_workspace = arcpy.GetParameterAsText(1)
arcpy.env.workspace = sde_workspace

try:

  #Converting Feature Classes to Shapefiles in Temporary workspace
        arcpy.FeatureClassToShapefile_conversion(["postgretest.postgres.boundaryline","postgretest.postgres.cartographicsymbol","postgretest.postgres.cartographictext", "postgretest.postgres.topographicarea", "postgretest.postgres.topographicline","postgretest.postgres.topographicpoint"], temp_workspace)

#Looping through each shapefile and updating the features in featuredataset in enterprise geodatabase
arcpy.env.workspace = temp_workspace
fcList = arcpy.ListFeatureClasses()

for fc in fcList:

    #Boundary Line............................................................................
    if(fc == "boundaryline.shp"):
        fc1 = fc
        fc2 = sde_workspace + "\postgretest.postgres.boundary_line"

        print "Starting for Boundary Line..."

        #set up cursor for adding features
        inscursor = arcpy.da.InsertCursor(fc2,["postgretest.postgres.boundary_line.fid","FEATURECODE","VERSION","versiondate","accuracyofposition","make","physicallevel","physicalpresence","style_description","style_code","SHAPE@"])
        with arcpy.da.SearchCursor(fc1,["fid_","featurecod","version","versiondat","accuracyof","make","physicalle","physicalpr","style_desc","style_code","SHAPE@"]) as innercursor:
            for searchrow in innercursor:
                geom = searchrow[10]
                row_values = [searchrow[0], int(searchrow[1]),int(searchrow[2]), searchrow[3],int(searchrow[4]),int(searchrow[5]),int(searchrow[6]),searchrow[7],searchrow[8], int(searchrow[9]),geom]
                inscursor.insertRow(row_values)

        del(innercursor)
        del(inscursor)
        arcpy.Delete_management(fc1)

Source Table:
"wkb_geometry" - Geometry
"fid" - Text Length = 0
"featurecode" - Long Numeric
"version" - Long Numeric
"versiondate" - Text Length = 0
"accuracyofposition" - Text
"make" - Text
"physicallevel" Numeric Long
"physicalpresence" - Text
"style_description" - Text
"style_code" - Nummeric Long


texts: Creating perpendicular line from points to existing line?

texts:
I would like to create perpendicular lines from points to line, similar to this:

I have ArcGIS 10.1, so "queryPointAndDistance" doesn´t work for me, and basic licence level, so I can´t use "Generate Near Table" neither.

I thought about following approach:

Coordinates X1Y1 are known (coordinates of the point feature)
acquire X2Y2 coordinates of the closest point on the line segment
create line feature from X1Y1 to X2Y2
new line feature is perpendicular to existing line

I tried to use "distanceTo" method from point feature to line feature, but I was unable to aquire coordinates of the closest point of the line segment (X2Y2), which should be used during "distanceTo" method.
Is it possible to get coordinates of the closest point on the line segment from "distanceTo" method or should I use another approach?

point_feature = "points_select"
line_feature = "line_select"

pointcursor = arcpy.da.SearchCursor(point_feature, ['SHAPE@'])

for row in pointcursor:
      geometry = pointcursor[0]
      linecursor = arcpy.da.SearchCursor(line_feature, ['SHAPE@'])
      for row in linecursor:
          newgeometry = linecursor[0]
          dist = newgeometry.distanceTo(geometry)
      print dist


texts: MOSAIC: Black grid lines between GEOTIFFs Using gdal_merge with Python

texts:
I've been trying to map an area of 35 images. Since they were in UTM I converted to Mercator projection using gdalwarp. I then tried to merge the images with the command:

gdal_merge -init 255 -o I:\newfilename.tif --optfile imagelist.txt

which produced a mosaic that had missing data where each image would have overlapped, similar to black grid lines.

I made sure that the order of images in my text file was arranged in order but can't figure out what else it could be. I'm completely new to coding so anything helps.


texts: Making a GeoJSON layer editable with the Leaflet.Editable plugin

texts:
I'm new to Leaflet and the Leaflet.Editable plugin, but I'd like to make GeoJSON layer editable. When I mean "editable", I mean a user can stretch the polygons at will.

Here is a good example of what I'd like.

My current code is as follows:

var map = L.map('mapid', {editable: true}).setView([39.08344, -77.6497145], 13);

L.tileLayer('Some Valid Tile Layer',
    {
    attribution: 'Myself',
    maxZoom: 18
}).addTo(map);

$.ajax({

    url: 'http://127.0.0.1:8000/static/map/geoJSON/elementary_boundaries.geojson',
    dataType: 'json',
    success: function(data){

        console.log('The data was received');
        L.geoJson(data, {

            onEachFeature: function(feature, layer){
                layer.enableEdit(); // incorrect, but code would go here

            }
        }).addTo(map);

    },

    failure: function(data){
        console.log('The data wasn\'t received');
    }

});

I don't fully understand features vs layers, and I fully don't understand leaflet's geometry construction flow, so I'm sure my problem is elementary.


texts: What's the best way to handle a 'Site Datum' in QGIS?

texts:
So I have a client with a site datum, which is to say they have an arbitrary coordinate reference system with coordinates between -1000 and 1000 in X and Y.

I have a list of ~20 coordinates that maps between a UTM projection and these coordinates, so I know I can calculate an affine transformation between them without too much difficulty.

My question is, can I handle this case with custom CRS, so that using the site datum is as simple as selecting it as the on-the-fly projection, or exporting data with the site datum transforms the data automatically?


texts: Oracle Spatial - SDO_GEOM.RELATE Inconsistencies

texts:
I'm using Oracle 12.2 (12c) and encountering some inconsistent results with the SDO_GEOM.RELATE function.
When I run the following code to check if two points are equal on the projected Irish grid, taking into account a tolerance, as expected I get a result of "EQUAL":
  SDO_GEOM.RELATE
  (
    MDSYS.SDO_GEOMETRY(2001,29902,MDSYS.SDO_POINT_TYPE(325511.001,313661.001,NULL),NULL,NULL),
    'EQUAL',
    MDSYS.SDO_GEOMETRY(2001,29902,MDSYS.SDO_POINT_TYPE(325511,313661,NULL),NULL,NULL),
    0.001
  )

But when I adjust to use a larger tolerance and a larger distance between points, I get a result of "FALSE":
  SDO_GEOM.RELATE
  (
    MDSYS.SDO_GEOMETRY(2001,29902,MDSYS.SDO_POINT_TYPE(325511.01,313661.01,NULL),NULL,NULL),
    'EQUAL',
    MDSYS.SDO_GEOMETRY(2001,29902,MDSYS.SDO_POINT_TYPE(325511,313661,NULL),NULL,NULL),
    0.01
  )

Here's a breakdown of my results for a given tolerance/distance:
0.0001 - FALSE
0.001 - EQUAL
0.01 - FALSE
0.1 - FALSE
1 - EQUAL

What is going on here?


texts:
I'm trying to do a gdal hillshade using this code:

  import os, subprocess
  #create hillshade using gdal
  subprocess.call(["gdaldem", "hillshade", r"/data/brent/SRTM DEM Training for Zamboanga Peninsula/raw_DEM.tif", \
     r"/data/brent/SRTM DEM Training for Zamboanga Peninsula/zamb_try_hillshade.tif", "-z 1.0", "-s 111120.0", \
     "-az 315.0", "-alt 45.0", "-of GTiff"])

It runs in the python console but there is no resulting tif file when I checked the ouput folder.

Did I miss something in my code? or there is a bug with gdaldem hillshade


texts: How to batch export using r.out.gdal command in GRASS

texts:
I created a list of files (NPPFILES) using g.mlist command in GRASS 7.
Now I would like to export these files as .tiff files.
I tried this command in GRASS prompt

for files in $NPPFILES; do
r.out.gdal input="$NPPFILES" output="$NPPFILES.tiff" format=GTiff ; done

This is not working for me. I'm running it in the directory where output files need to be stored, still no result.


texts: Converting Raster Density File into Point using ArcGIS Desktop?

texts:
I used the Point Density tool to create raster files. I plan on scaling these density values by dividing them by the greatest density value and then multiplying the result with a certain value that is associated with that feature class. But first I am planning on converting each raster cell(which holds a density value) into a point. This way, I can manipulate the values in the attribute table.

I am trying to use the Raster to Point tool but keep running into this issue?

I'm not quite sure what it means when ArcMap says the output point features are empty.


texts: What this error means on Google Earth Engine console? "Image (Error) Invalid JSON"

texts:
When I run my code on GEE I receive the following error in the line print(unitScale).
Image (Error) Invalid JSON.
var features = (elevation
    .addBands(slope)
    .addBands(curve)
    .addBands(landcover) // Can disaggregate to dummy features in Python
    .addBands(twi)
    .addBands(spi)
    .addBands(sti)
    .addBands(mtpi)
    .addBands(lith)
    .addBands(GSWPerm)
    .addBands(GSWDistSeasonal)
    .addBands(aspect)
    .addBands(hand)
    .clip(region))
// .clip(elevation.geometry()));

features = features.reproject({
    crs: aspect.projection().crs(),
    scale: 30
}).float()

var minMax = features.reduceRegion({
    reducer: ee.Reducer.minMax(),
    scale: 30,
    maxPixels: 40e9,
    bestEffort: true
})

// use unit scale to normalize the pixel values

var unitScale = ee.ImageCollection.fromImages(
    features.bandNames().map(function (name) {
        name = ee.String(name);

        var band = features.select(name);

        return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))))
    })).toBands().rename(features.bandNames());

print(unitScale) /// I receive the error in this line ///

print(features)


texts: Using python Point-In-Polygon operation to join grid with points based on points falling within grid?

texts:
I want to know how to select values from the xarray DataArray based on the location (geo_df.geometry) and time (geo_df.plant_date & geo_df.cut_date) of rows in the geopandas GeoDataFrame. I want to join them as 'features' in an output GeoDataFrame.

My datasets:

Packages I'm using:

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely import geometry
import xarray as xr

I have a geodataframe storing lat/lon POINTS which corresponds to households. The index column is the id of the households.

geo_df.head()

Out[]:
  crop_name     xxx     cut_date plant_date                       geometry
0   SORGHUM  0.061029 2011-11-10 2011-11-10 POINT (37.89087631 14.35381619)
1    MILLET -0.104342 2011-10-19 2011-10-19 POINT (37.89087631 14.35381619)
2   SORGHUM -0.031697 2013-11-26 2013-11-26 POINT (37.89087631 14.35381619)

I have an xarray object storing GRIDDED vegetation health data (NDVI).

ndvi_df = xr.open_dataset(geo_data_dir+ndvi_dir).ndvi

Out[]: <xarray.DataArray 'ndvi' (time: 212, lat: 200, lon: 220)>
[9328000 values with dtype=float32]
Coordinates:
  * lon      (lon) float32 35.024994 35.074997 35.125 35.174988 35.22499 ...
  * lat      (lat) float32 14.974998 14.924995 14.875 14.824997 14.775002 ...
  * time     (time) datetime64[ns] 2000-02-14 2000-03-16 2000-04-15 ...
Attributes:
    long_name:   Normalized Difference Vegetation Index
    units:       1
    _fillvalue:  -3000

I have a geodataframe storing a POLYGON which corresponds to a country.

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
ethiopia = world.loc[world["name"] == "Ethiopia"]

Visual Summary:

My datasets plotted on top of one another look as follows (plotted annually for demonstration purposes).

(ndvi_df.loc[f'{year}-01-16T00:00:00.000000000':f'{year}-12-16T00:00:00.000000000']
 .mean(dim='time')
 .plot(cmap='gist_earth_r', vmin=-0.1, vmax=1)
)

ax = plt.gca()

ethiopia.plot(alpha=0.2, color='black', ax=ax)

(geo_df
 .loc[ (lsms_geo_1["cut_date"] > f'{year}-01-01') & (lsms_geo_1["cut_date"] < f'{year+1}-01-01') ]
 .plot(markersize=6 ,ax=ax, color="#FEF731")
)
ax.set_title(f'{year} Mean NDVI and Households')
plt.show()

Ideal Output:

I want as an output, a geodataframe with extra columns telling me the NDVI values in the PRECEDING MONTHS for the pixel which the households are inside.

The index column is the id of the households.

like this:

  crop_name     xxx     cut_date plant_date                       geometry  ndvi_month_0  ndvi_month_1  ndvi_month_2
0   SORGHUM  0.061029 2011-11-10 2011-11-10 POINT (37.89087631 14.35381619)          0.3           0.3           0.3
1    MILLET -0.104342 2011-10-19 2011-10-19 POINT (37.89087631 14.35381619)          0.6           0.6           0.6
2   SORGHUM -0.031697 2013-11-26 2013-11-26 POINT (37.89087631 14.35381619)          0.1           0.1           0.1

I would also like to know how to subset my data in xarray object by using the geodataframe polygon ethiopia.

NOTE: Reposted from stack overflow here because it seems a GIS related question.


texts: How to use pythons' GeoAlchemy2 to filter on dynamically generated Points?

texts:
I'm trying to do something like this SQL query in Geoalchemy2:

SELECT district_polygons.district
FROM markers, district_polygons
WHERE ST_Contains(district_polygons.multi_poly, Point(markers.marker_long, markers.marker_lat))

This query works great as is in MySQL. This shows me which district each marker is in.

I want to replicate something like this using Geoalchemy2 / SQLAlchemy. This is the what I tried:

results = dbSession.query(District_Polygon.district).\
    filter(func.ST_Contains(District_Polygon.multi_poly, WKTElement('POINT({0} {1})'.format(Marker.marker_long, Marker.marker_lat)))).all()

Unfortunately the formatter isn't working, and sqlalchemy output looks like this:

SELECT district_polygons.district AS district_polygons.district
FROM markers, district_polygons
WHERE ST_Contains(district_polygons.multi_poly, ST_GeomFromText(%s))

INFO sqlalchemy.engine.base.Engine ('POINT(markers.marker_long markers.marker_lat)')
[]

How can I replicate the functionality of the SQL command in the ORM model of Geoalchemy / SQLAlchemy? I'd like to avoid raw SQL if I can.


texts: OpenLayers code didn't response popup featureinfo

texts:
Hey I'm trying to popup the featureinfo of WMS layer using OpenLayers code. But it does'nt work. Instead the cursor moves in the state of continous moving when I click on map. Here is my code:

<title>Untitled 1</title>
        <script src="http://localhost:8080/geoserver/openlayers/OpenLayers.js" type="text/javascript">
    </script>
    <script type="text/javascript">
        var map;
         var wms;
         var wms1;

        function init(){

            alert("inside function");
            var bounds = new OpenLayers.Bounds( 143.853, -43.648,
                148.479, -39.574);

        var options = {
        controls: [],
        maxExtent: bounds,
        maxResolution: 0.0180703125,
        projection: "EPSG:4326",
        units: 'm'
        };

            map = new OpenLayers.Map('map', options);
            wms = new OpenLayers.Layer.WMS("Peoples Colony",
                               "http://localhost:8080/geoserver/wms",
                               {
                                   LAYERS: "topp:tasmania_state_boundaries",
                                   STYLES: '',
                                   format: "image/jpeg",
                                   transparent: "true",
                                   tiled: 'true'

                               },
                               {

                    isBaseLayer: true,
                    yx : {'EPSG:4326' : true}

                               }
                               );
                               wms1 = new OpenLayers.Layer.WMS("tasmania_roads",
                               "http://localhost:8080/geoserver/wms",
                               {
                                   LAYERS: "topp:tasmania_roads",
                                   STYLES: '',
                                   format: "image/jpeg",
                                   transparent: "true",
                                   tiled: 'true'

                               },
                               {

                    isBaseLayer: false,
                    yx : {'EPSG:4326' : true}

                               }
                               );
        map.addLayers([wms,wms1]);
        map.zoomToMaxExtent();

             map.addControl(new OpenLayers.Control.Navigation());
              map.addControl(new OpenLayers.Control.Scale($('location')));

        var info = new OpenLayers.Control.WMSGetFeatureInfo({
        url: "http://localhost:8080/geoserver/wms",
        title: 'Identify features by clicking',
        queryVisible: true,
        eventListeners: {
            "getfeatureinfo": function(event) {
                map.addPopup(new OpenLayers.Popup.FramedCloud(
                    "chicken",
                    map.getLonLatFromPixel(event.xy),
                    null,
                    event.text,
                    null,
                    true
                ));
            }
        }
    });
    map.addPopup(info);
    info.activate();

}

</script>

Anyone plz help me out? I think that I didn't mention the layers here. Is this the problem?


texts: Python syntax error with If/elif statement in field calculator when using concatenation

texts:
I am unsure what syntax issue I am having with my If/Elif statement. I want to not use a column in my calculation if the value is blank. The end goal is to not have extra spaces before, after and between values.

Please be kind. I am still learning Python for scripting in ArcGIS Desktop.

Here is my Prelogic and Expression.

Prelogic:

def ifFields(SubAdd,LandM, ZipN):
  if SubAdd== " ":
    return "!LandmarkName! + ' ' + '('+ !ZipName! + ')'"
  elif LandM== " ":
    return "!SubAddress! + ' ' + '('+ !ZipName!+')'"
  elif ZipN== " ":
    return "!SubAddress! + ' ' + !LandmarkName!"

Expression:

ifField(!SubAddress!, !LandmarkName!, !ZipName!)


texts: How can I determine which service is associated with a particular ArcSOC.exe process and vice versa?

texts:
I've got a development arcgis server machine configured with multiple mapservices, gp services and SOE's.  Sometimes an arcgissoc.exe process gets hung.  There also seem to be too many SOC's running based on what I can see in ArcCatalog.

How can I figure out which ArcSOC.exe (process id) goes with which service in arccatalog?

I'm running 10.0.


texts:
When I open a HDF file in ArcMap, the display of HDF file is shown as Fig.1. The HDF file has many bands, Arcmap will popups the window when open the HDF file regularly.

However, there is no pop-up window  when I open HDF file in the Arcmap software and the display is shown as Fig.1. I don't know the reason for the error. Have anyone encountered the error before?

                                          Fig.1

Now,I think it a lot,maybe the popup window failed to be opened,so the ArcMap just opened the first band of HDF file.If it is true, how I make the popup window be opened?


texts: How to perform trilateration using 3 lat/lon points without distances?

texts:
I read this thread Trilateration using 3 latitude and longitude points, and 3 distances explaining about trilateration using 3 latitude and longitude with measured distance.

The problem is I want to find out an unknown target location with only 3 knowns latitude and longitude co-ordinates without knowing each distance point. For example, I have 3 following longitudes and latitudes:

  lat="-6.28530175" lon="106.9004975375"
  lat="-6.28955287" lon="106.89573839"
  lat="-6.28388865789474" lon="106.908087643421"

Thank you very much :)

Based on @D.E.Wright answer I am googling to find out more and I get this following forum http://www.coderanch.com/t/453432/Programming/implement-cell-triangulation-mobile-phones, it gives 3 steps to triangulation. The problem I couldn't understand what he means about "Calculate distance from the first tower based on speed which gives a radius value".  And with @Kirk Kuykendall suggestion I decide to put my 3 lat/long to google map and I get this following image :

Where :
A, B, C : Cell Towers location

A : lat -6.28955287, lon 106.89573839 with cell tower Radius 6000m
B : lat -6.28530175, lon 106.9004975375 with cell tower Radius 6000m, and
C : lat -6.28388865789474, lon 106.908087643421 with cell tower Radius 6000m

Now, with these data can I doing a triangulation to get my real position? And how is the mathematical method? I've try to modify phyton code to remove the distance variable from @wwnick in this thread Trilateration using 3 latitude and longitude points, and 3 distances but until now I couldn't find the answer.


texts: If/Then Statement in ArcMap to Remove Duplicates

texts:
I need help writing an if/then statement in ArcMap 10.5.

Scenario:
I generated a near table, but since the closest features (polylines) were mostly all overlapping with the input (multipoints), the table "ranked" the outputs as tying since the distances were zero.

I would like to re-rank these values so that the rank # is always unique, as I need to convert this near table into a pivot table next. It doesn't matter what ranking they get, as they're technically all the same distance.

Any ideas for an expression to rewrite the ranking # if there are duplicates, or a different workaround?

Example:
IN_FID 2 found three features that were within my search radius. Instead of listing them as rank 1, 2, 3 it's listing it as 1, 1, 2 because the Near_Dist is 0. I just want a way to re-rank these so that the table doesn't have rank numbers that tie, like how this example has two "1"s. These can be arbitrary numbers, and don't need to be in this field. I just need a field to give a unique number between 1 and 3 for every IN_FID in this Near Table.


texts: How to automatically merge Texas survey block data

texts:
I have Data from the Texas rail road commission.  I'm trying to get the block data to display with a thick line.  The problem is the data is in surveys.  I exported a copy of the shapefile.  If you take a look at the L2BLOCK colunm you see that I started merging them together, but it take forever if you have to do so for a lot of counties.

I want the blocks to have a thicker line to seperate them from the surveys.  Here's what the table above looks like.

You can see that some of them are merged and some not.  Is there way to automate merges based on a same block numbers in the attribute table?


texts: I'm using QGIS 2.18.4, I wanted to add "scale text" but can't find the tool for adding the scale text on the composer. I only see an option of add text and not "scale text" as in other GIS Software like ArcGIS. Is there anyone who knows this?


texts: Wrong result using virtual layer sum with join and group by

texts:
I am trying to join the grouped sum of two layers using virutal layers in QGIS 3.2.

My query looks like:

select nvewi.zone as "Zone-Type", sum(nvewi.bev) as "Inhabitants Isochrones", sum(nvewp.bev) as "Inhabitants Puffer"
from nvewi
inner join nvewp
on nvewi.zone=nvewp.zone
group by nvewi.zone

Table nvewi:

bev zone
50  a
20  a
30  a
10  b
20  b
30  b
2   c
2   c

Table nvewp:

bev zone
10  a
10  a
20  a
5   b
3   b
2   b
1   c
1   c

The result I get is obviously wrong:

Zone-Type   Inhabitants Isochrones  Inhabitants Puffer
a   300 120
b   180 30
c   8   4

The result how it should be:

Zone-Type   Inhabitants Isochrones  Inhabitants Puffer
a   100 40
b   60  10
c   4   2

Where is my mistake? And how to perform this join correctly getting the correct sums?


texts: Reproducing Difference operation from QGIS in pure Python?

texts:
I have two GeoJSON files and I want to make the difference of those two layers, QGIS does it really well, but I need to do it in pure Python 2.7.

Anyone know how to do it ?


texts: Is it necessary to use projected coordinate system for sf (function: st_nearest_feature) package?

texts:
I have two csv files, each with lat/long coordinates. I converted them into sf objects, giving each GCS (geographic coordinate system) crs. When I use the st_nearest_feature function, I get this warning:

although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar

I'm just getting the index of the closest points and don't need the exact distance.

Are there any problems in not converting the coordinate systems ahead of time that I should be aware of?


texts: How to use map units to define radius of point feature in styles in QGIS?

texts:
I have a point feature in a map which is in EPSG32723 (WGS84 - UTM zone 23s) and I need to style it as follows: A center point in red and a circle below it which has a particular size in meters.

I tried doing that in styles using map units for the size of the second but when I measured it using measure tool the circle had half (3 meters) of what it was supposed to be (6 meters).


texts: Alternative to Google Earth to edit and organize KML files

texts:
I like the way Google Earth displays KML as a "scene graph", that is, a tree structure in the left panel.

This allows you to rearrange features in folders, change layer order, etc. via drag-and-drop.

However, this (rather useful) functionality is a bit too basic in Google Earth. For example, google "imports" the KML in its own structure, and if you want to "save back" the changes, you have to manually point to the original file in order to overwrite it.

Also, upon save, Google Earth just overwrites the styles contained in the original file (pretty much what Inkscape does with SVG). That's pretty rude :o(

So, my question is: are there other alternatives, specifically using the Tree View representation with drag-and-drop rearrangement? I don't even need map display (although that would be very useful)


texts: Creating null values when creating raster in arcgis

texts: I am converting a large set of polygons into a single raster where all of the polygon values are 1. When I do this I do not have null values (0) for the rest of the cells.  I looked through polygon to raster and can't find a set null value.  I tried setnull within raster calculator and it just changes the values of the polygons.  What I need is a raster with 1 for protected area and 0 for no protected area so I can run a number of calculations.


texts: QGIS Composer - Read entered textbox value in another textbox

texts:
I'm Limited to QGIS 2.18 composer for several reasons.
I have two textboxes.
BOX1 - User INPUT value, where user enters an integer value.
BOX2 needs to read this value to be able to calculate results.
Can I use the item_uuid to read value in BOX1 to calculate a result in BOX2?
If so, what expression would this be?


texts: gdal.Warp cutline creates empty file when run on the server, but works perfectly fine through Python command line

texts:
Trying to cut geotif with a shp file.

When I open Python command line interpreter and run

from osgeo import gdal

gdal.Warp(outfile, infile,format='GTiff', cutlineDSName=shpfile,  srcSRS='EPSG:4326', dstSRS='EPSG:4326', dstNodata=-9999)

it creates the geotif as expected, though shows following output:

Warning 1: for band 1, destination nodata value has been clamped to 0, the original value being out of range.
<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x7f8c72ac5b10> >

When I do the same within a web server, doing file preparation (projection, transparency, etc.) first with gdal and then trying to cut the file with the same command I get an empty file (0 bytes) and server (gunicorn) restarts.

Machine:

Linux Mint 18.3

Libraries versions:

gdal==2.2.2
proj==4.9.3
geos==3.6.2
python==2.7.12
gunicorn==19.7.1

Does anyone have any idea how to make it working in the web server environment?


texts: Showing directional line symbol in ArcGIS for Desktop?

texts:
Is this possible in ArcGIS to have variable number of arrows on the line to show direction?

This is an ArcGIS example, where I can have fixed (1,2,3,4,5...) number of arrows and it stays the same doesn't matter what the line length is. First and last always stick to the end, doesn't help with network mapping, where there is usually the node to show. Having just one in the middle is never enough, check line at the right...

No such issues with ArcView 3:


texts: Generating DEMs from point clouds and differencing them?

texts:
I have two point clouds obtained via structure from motion (ground photos of sand-dunes + Visual SfM/Agisoft Photoscan).

If I wanted to generate a DEM from these point clouds and difference them to look at elevation/sediment change, how would I go about doing this?

Is there any free software that can achieve this easily? I am a beginner at using GIS software.

At the moment I am using the free trial of Agisoft Photoscan, which can export the point cloud in Wavefront OBJ, Stanford PLY, XYZ text file, U3D, ASPRS LAS, and PDF formats.


texts:
I would like to program a standalone application using C++. In version 3.4 LTR, my program ran fine and I did not receive any notification when running it. However, when I compiled using the qgis library version 3.10.3 LTR, I received the following warnings:

proj_get_authorities_from_database: Cannot find proj.db
proj_create_from_wkt: Cannot find proj.db
proj_create_from_database: Cannot find proj.db
proj_create_from_database: Cannot find proj.db
proj_create_from_database: Cannot find proj.db

Can someone help me?


texts:
Is it possible to export a QGIS project to Mapnik xml?

I'm using QGIS 2.4 and Quantumnik does not seem to work with that version.


texts: Add a 'cloud' band to Sentinel 2 images (GEE)

texts:
This script is used to mask out cloud in sentinel 2 images using a threshold for band 4:

// Sentinel 2 cloud function
var cloudfunction_ST2 = function(image) {
  // If band 4 is higher than 2500, the pixel is considered as cloudy.
  var b4 = image.select("B4");
  // Get pixels above the threshold.
  var cloud = b4.gt(2500);
  // Create a mask from high likelihood pixels.
  var cloudmask = image.mask().and(cloud.not());
  // Mask those pixels from the image.
  return image.updateMask(cloudmask);
};

My question is, how to add a 'cloud' band to the image collection?


texts:
I want to get the data within BBOX and form the CQL query URL as following:

http://[localhost]:8080/geoserver/mymap/ows?service=WFS&version=1.0.0&request=GetFeatur
e&typeName=layer:layername&outputFormat=json&BBOX=-70.95699,50.89799,-71.17946,52.04493,EPSG:4326&cql_filter=(name+like+%27%25h%25%27%29


texts:
I'm trying to upscale/downsample some elevation raster data at 2m resolution to 20m to speed up my analysis. To do this I'm reading in the individual files, assigning their crs and upscaling them using the 'out_shape' parameter in rasterio's read function.

The issue I'm encountering is that whilst this reduces the resolution of the image it doesn't appear to be doing so whilst accounting for the crs, this means that my tiles then have huge gaps which match the ratio I'm upscaling with. The unscaled image can be seen on the left and the scaled one on the right.

Code:

import geopandas as gpd
import rasterio as rio
from rasterio import features

results = []
test_files = ['tl9001_dtm_2m.asc', 'tl9002_dtm_2m.asc']

for test_file in test_files:
    with rio.open(f'data/test/{test_file}', crs='epsg:27700') as src:
        src_meta = src.meta.copy()
        src_affine = src_meta.get('transform')

        band = src.read(1, masked=True, out_shape=(50, 50)) # Input is 500 x 500

        for geometry, raster_value in features.shapes(band, transform=src_affine):
            result = {'properties': {'raster_value': raster_value}, 'geometry': geometry}
            results.append(result)

gpd_results = gpd.GeoDataFrame.from_features(results)
gpd_results.plot(column='raster_value')

Data: https://drive.google.com/drive/folders/18fyq8TpZgRvScDWbsvzn8axtKJkFacAE?usp=sharing