texts: clipping a vector layer with another vector layer using python

texts:
Sorry if this there is a really obvious solution to this question. I did not have any luck googling and looking at different books.

I have an OSM vector layer that I want to clip with another vector layer for the boundary of a city. Now I can do this easy enough in QGIS, but I wanted to do this in python since I have a collection of images to process and clip.

I imagined that clipping is such a common operation, that there must be clip functions in the major geospatial packages, but I have not been able to find anything like that so far. I was kinda surprised, or was wondering if I was looking under the wrong headings, terminology?

There are a few posts on this topic listed below.

Clipping raster with vector layer using GDAL

Clipping shapefile with Python?

"Clip" hole into Shapefile

Now the first post is about 7 years old, and the later ones are a couple of years old.

Is the best way to do this to use something like geopandas and to then keep all geometries that intersect with the city polygon? But then how do I crop OSM polygons that are not completely contained within the city polygon boundary. I would like to keep partial polygons.

Is the code in the links above still valid, or are there better ways to do this now?


texts: Create grid with a prespecified number of cells inside a shapefile with R

texts:
I need to create a grid inside a shapefile with a fixed the number of cells n (inside the shapefile). I tried a few stuff see bellow but nothing worked so far. Any ideas?

The data shapefile can be found here:

https://geodacenter.github.io/data-and-lab/scotlip/

library(maptools)
library(sp)
library(raster)
library(rgeos)

scotlip <- readShapeSpatial("~/Desktop/scotlip/scotdistricts.shp")
plot(scotlip)

n <- 100  # the number of cells I want include the shapefile.

# Attempt 1

set.seed(12)
sppoint <- spsample(scotlip, n = n, type = "regular")
nrow(sppoint@coords)
# first problem: 92 instead of 100 points

points(sppoint, col = "red", pch = 19, cex = .5)

grid <- SpatialPixels(sppoint)
grid <- as(grid, "SpatialPoygon")
plot(grid)
plot(scotlip, border = "red", add = T)
# second problem: it doesnt cover the shapefile, and thus the areas differ

# Attempt 2

xmin <- scotlip@bbox[1,1] - 10000
xmax <- scotlip@bbox[1,2] + 10000
ymin <- scotlip@bbox[2,1] - 10000
ymax <- scotlip@bbox[2,2] + 10000

ext <- extent(c(xmin, xmax, ymin, ymax))

# try fixing the resolution based on the number of grid cells required to cover the area, but
# this doesn't work too, because of the grid cells that intersect with the land.

tot.area <- sum(area(scotlip))

res(r) <- sqrt(tot.area/n)
grid <- SpatialPixelsDataFrame(SpatialPoints(coordinates(r)),
                               data = data.frame(layer = 1:nrow(coordinates(r))))

grid <- as(grid, "SpatialPolygonsDataFrame")
plot(grid)
plot(scotlip, border = "red", add = T)

int <- intersect(grid, sctolip)
sum(area(int)) - sum(area(scotlip))
plot(int)
nrow(int)


texts: Data loss after compressing geodatabase in ArcSDE

texts:
I suffer sometimes from data loss after compressing ArcSDE geodatabase, this  geodatabase contains versioned feature classes, so the compressing script reconciles and compress the geodatabase:

# Execute the ReconcileVersions tool.
arcpy.ReconcileVersions_management(db_path, "ALL_VERSIONS", "sde.DEFAULT", versionList, "LOCK_ACQUIRED", "NO_ABORT", "BY_OBJECT", "FAVOR_TARGET_VERSION", "POST", "DELETE_VERSION", "c:/reconcilelog.txt")

# Run the compress tool.
arcpy.Compress_management(db_path)

it seems that the main reason is a bad use of ReconcileVersions_management, what are the parameters recommanded ? and how i can test if data loss would happen again ?

the geodatabase is based on ArcSDE 10.1 and Postgre SQL 9.0.5.


texts: Is there an open source tool or module in python that does parcel division of polygons like in ArcEditor?


texts: How to make the coordinates of point appears as labels?

texts:
I am working on a school project that requires us to have a map with the points indicating different station with its coordinates appearing as a label.  I just don't know how to do it?

I think it could be possible but we just started using QGIS so I know very little.


texts: Dividing polygon into tiny polygons of particular size

texts:
I would like to divide/cut an irregular polygon into tiny polygons of a particular size(1.6m x 1m) in such a way that most of the irregular polygon area has to be utilized (an OPTIMIZATION MODEL)
The length and width of the polygon can be interchanged (either 1.6m X 1m (or) 1m X 1.6m)
So, in the end, I need to have as many polygons of size (1.6m X 1m) as possible.

I used the "Create Grid" feature but it just cuts the whole polygon in a particular fashion.

But what I also want is that here, a blue polygon can also be cut in a vertical manner(1m x 1.6m) too.

So, I would like to know where there is a plugin for this in QGIS/ArcGIS or any python script for this kind of polygon optimization?


texts:
There's loads of web-based tile layers out there which are great in many ways, but one problem I'm having is making them look nice when exporting a map to print or as a jpeg.

The problem consistently seems to be that they display too much detail, which is fine for working with the maps on-screen, but for most publications I want to simplify the map and have everything nice and clear (particularly thinking of boundary and place-name layers). It seems to be that one way to achieve this would be to overrule the tile layer's Levels of Detail/Scale settings to get a lower Level of Detail despite the map being a large scale.

Is that possible (for a third-party, web-based tile layer)?

If not (or even if so) does anyone have an better ideas for achieving a similar effect?

I use both QGIS and ArcGIS Pro so a workaround in either would be OK.


texts:
I am a novice in arcgis.

Initially I had a raster of land use with several classes which I used a neural network model to predict future land use based on one class which is agriculture, in the study area. Processing the raster using arcmap and the neural network finally gave me a raster with the predicted future land use scenario with two classes (the one being modelled =1, and the remaining land use which the model classified as 0). Now I want to revert back to a raster with the original classes (of course including the agriculture that has changed).

These are images of the two rasters.

[2


texts:
I am trying to create symbol in QGIS using svg. When I place the point in QGIS it symbology remains same.

Is there any possibility to stretch the symbol in QGIS?
Can we make stretchable point symbols in QGIS?


texts: Joining .csv to attribute table in ArcMap?

texts:
Is there an ArcMap 10.5 method of joining an excel file (.csv) to a shapefile/attribute table using several fields instead of an explicit identifier?

I've been having trouble with this, and joins simply appends it and fills all rows in the attribute table with the first row from excel.

They don't share a specific number ID, but both have columns (fields?) such as "Township", "Lot", and "Concession" that share the same data. I am looking to join it using these as identifiers and then fill the remaining columns with the excel data.

In summary, is there a way to simply fill the fields in the attribute table with the data from excel?


texts: create buffer around locator-derived point

texts:
I am trying to create a buffer around a point derived from a locator search, using this example:
http://help.arcgis.com/en/webapi/javascript/arcgis/help/jssamples_start.htm
..but I am getting no results.  I am not sure how to pass the correct information from the locator result graphic to the buffer functions (sorry if this has been answered before, i couldn't find a solution).

Thanks, code is below..Jason

    ////////////// locator func ///////////////////////////////////

  function showLocatorResults(candidates)
  {
    var candidate;
    var symbol = new esri.symbol.SimpleMarkerSymbol();
    symbol.setStyle(esri.symbol.SimpleMarkerSymbol.STYLE_SQUARE);
    symbol.setColor(new dojo.Color([153, 0, 51, 0.75]));
    var geom, graphic;

    dojo.every(candidates, function (candidate) {

        if (candidate.score > 80)
        {
           // extent code removed for brevity
            var attributes = { address: candidate.address, score: candidate.score,    locatorName: candidate.attributes.Loc_name };
            geom = candidate.location;
            graphic = new esri.Graphic(geom, symbol, attributes);

            map.graphics.add(graphic);
            return false;
        }
    });

    doBuffer(graphic);
    }

     ////////// buffer funcs ////////////////////////////////////////////////
     function doBuffer(g) {
    //map.graphics.clear();

    bufferParams = new esri.tasks.BufferParameters();
    bufferParams.distances = [100];
    bufferParams.geometries =  [g.geometry] ; // problem here?

    bufferParams.unit = esri.tasks.BufferParameters.UNIT_FOOT;
    bufferParams.outSpatialReference = map.spatialReference;
    geometrySvc.buffer (bufferParams, showBuffer);
}

function showBuffer(geometries) {

  var symbol = new esri.symbol.SimpleFillSymbol(
    esri.symbol.SimpleFillSymbol.STYLE_SOLID,
    new esri.symbol.SimpleLineSymbol(
      esri.symbol.SimpleLineSymbol.STYLE_SOLID,
      new dojo.Color([0,0,255,0.65]), 2
    ),
    new dojo.Color([0,0,255,0.35])
  );

  dojo.forEach(geometries, function(geometry) {
    var graphic = new esri.Graphic(geometry,symbol);
    map.graphics.add(graphic);
  });

}


texts: Executing arcpy.createTable_management with parameters from arcpy.GetParameterAsText?

texts:
I'm trying to create a script tool in ArcMap that creates an attachment match table.  The script fails giving this:

Traceback (most recent call last):
  File "H:\py\PhotoMatch_test.py", line 40, in <module>
arcpy.CreateTable_management(outputGdb, matchTable)
  File "c:\program files\arcgis\desktop10.2\arcpy\arcpy\management.py", line 14419, in CreateTable
raise e
ExecuteError: ERROR 999999: Error executing function.
Failed to execute (CreateTable).

My code looks like this:

import arcpy, sys

inputFeatures = arcpy.GetParameterAsText(0)
picfolder = arcpy.GetParameterAsText(1)
matchTable = arcpy.GetParameterAsText(2)

outputGdb = arcpy.env.scratchGDB
arcpy.CreateTable_management(outputGdb, matchTable)

The parameter properties of the script tool:

I'm obviously a beginner.


texts: Determining which "color" each image band represents?

texts: I have a satellite image that has a few bands in it. But I don't know what colors they are supposed to represent, nor do I know where to look for this information--the file's metadata say nothing about this. Is there a way to determine what color each band represents? How can I find out which band is red, blue, NIR, etc.?


texts: How to update country borders in matplotlib.basemap?

texts:
I'm using basemap-1.0.7 and generated an image of Africa with the drawcountries() method.

The border between Sudan and South Sudan is missing. Does someone know how to update this country borders?


texts: Split Multiple Polygons into equal squares?

texts:
I have 1 single shapefile of about 70k square polygons in a grid.  I want to run a process/tool/script to split all 70k polygons into 4 equal parts.  The polygons are contiguous, so each polygon only has 4 vertices.  Basically I have a shapefile of sections for a sate, and I want to split it into 1/4 sections so each original square section is broken into 4 equal 1/4 parts.

How can I split all 70,000 contiguous square polygons at once with only 4 vertices each into 4 equal parts to create 4 polygons within each original polygon?


texts:
I've generated a Geospatial PDF with GDAL 1.11.3 library.

When I open it with Global Mapper v18, it looks upside down and it's not well georeferenced.

I've opened it with QGIS and with Acrobat Reader, and it seems it's ok.

So I don't know if it's a bug, my file is wrong or Geospatial PDF is not supported.

This is my file geomapa.pdf


texts:
I am running the following line in a bash script:
gdal_translate -of JP2KAK ${2}.vrt ${2}.jp2
to convert a virtual mosaic created with gdalbuildvrt to a JP2 file of the same name. The command produces the output with the artefact shown in the attached image. I have also tried using different QUALITY=n settings with the same result.  Creating a GeoTIFF from the same virtual mosaic works fine with no artefact.
Can anyone explain why the JP2KAK driver is causing the artefact shown in the attached image?


texts:
How to generate mbtiles in geoserver? using openlayers for displaying geoserver layer

for example calling wms layer like this

       new OpenLayers.Layer.WMS("Kanpur", "http://localhost:8080/geoserver/wms",
        {
                LAYERS: 'sample_data_old:sample',
                STYLES: '',
                format: 'image/jpeg',
                tiled: false,
                transparent: true,
                tilesOrigin : map.maxExtent.left + ',' + map.maxExtent.bottom,
                visibility:false
                },
        {
                isBaseLayer: true
        }
);

how can i call wps layer for creating mbtiles?

please refer this link for my question

http://docs.geoserver.org/stable/en/user/community/mbtiles/output.html


texts: Repeat only part of model with a while loop

texts:
The first part of my model isolates roads in a jurisdiction and then the second part of my model selects the roads in the jurisdiction that are greater than 500m and splits them.

I've input a while iterator to continue if any of the roads are greater than 500m but I'd like it to only apply to the second part of the model and ignore the first part. Is this possible?

EDIT: In order to simplify things I split the section I'd like to include in the while loop into its own section. The problem now is the input is not being updated so after each iteration the number of sections over 500m resets.


texts:
I guess mine is a very simple task, however I don't know how to solve it.
I m trying to display the outlines of a raster created from a polygon shp.
What I would like to come up with is an image where I can visualize the outlines of each poligon within my raster image. This image has two different colours (namely two "classes") and one is switched off, so only one is visible. Here is an example of what I need (but this is a shp, my image now only displays the yellow parts, while I also need the black lines):

I work both with ArGIS 10.1 and QGIS 1.8. Any idea?


texts: Generating a Placemark (Google Earth) via information link

texts:
I'm using the ArcGIS JavaScript API.

I'd like to have an infowindow pop up that has a link to a downloadable KML placemark. Is this possible without pre-generating thousands of individual KML files?


texts:
I have a gpkg that is populated exclusively by .jpeg files. I am attempting to
use the gdaladdo to add overviews into that gpkg. However, gdaladdo adds some of them as .jpeg, and adds the rest as .png. I was wondering if there is any way to configure gdaladdo so that it will only add the overviews as
.jpeg files.

The commands I run are:

gdal_translate -of GPKG myTiff.tif myGpkg.gpkg -co
TILING_SCHEME=InspireCRS84Quad -co ZOOM_LEVEL_STRATEGY=LOWER -co
TILE_FORMAT=JPEG

gdaladdo myGpkg.gpkg 2 4 8 16 32 64 128


texts: Using a list to access database from a layer

texts:
I am trying to access the database using a list in order to put it in a loop but I can't make it work. I want to replace the months with a list variable.

This is the hard core example:

   arcpy.GACreateGeostatisticalLayer_ga(in_ga_model_source="Areal
Interpolation",in_datasets="London X=London.Shape Y=London.Shape
F1=London.february",out_layer="february_interpolation")

I have list with the months which is:

    myMonths =['january', 'february',  'march','april','may','june','july',
'august', 'september', 'october',   'november', 'december']

I am trying to find a way to change the months within a loop. I tried the following but I am sure that it doesn't work.

    arcpy.GACreateGeostatisticalLayer_ga(in_ga_model_source="Areal
Interpolation",in_datasets="London X=London.Shape Y=London.Shape F1=London.%d"
 %myMonths[1],out_layer=myMonths[1]+"_interpolation")

I tried to find the answer but I was unable to find suitable search terms.


texts: Merge Multiple Tables into a New Table in PostGIS

texts: I am looking to merge a number of individual tables into a new table in PostGIS. This is an easy task when working with Shapefile data, but I am unsure how to do this in PostGIS.  Any help would be greatly appreciated.  I think I use Append, but I am unsure even where to start.


texts: Extracting DEM from imagery using ArcGIS for Desktop?

texts:
I'm using ArcGIS for Desktop and I'd like to know how can I extract the elevation model from a imagery of photos.
I have the have each image (jpg) geo referenced. After, I don't know what tool may use to analyst these pictures and get an 3d model or the DTM or DEM information ( something like Photomodeler Scanner does, a Dense Surface).

Can I do this with ArcGIS for Desktop (ArcMap)?


texts: Generating a SQL statement from a text file to make a selection

texts:
I have a text file of project names used in various workflows, it is structured as the following:
Project_2020_1m_20200212_0900
Project_2020_1m_20200213_0900

I want to use the text file to create an SQL statement, and then use that statement in a SelectLayerbyAttribute_management, and then CalculateField_management. I can create the statement but I don't seem to be able to make the selection (calculate fields is occurring on all the projects).
#create the sql selection string
datasetslist = datasets.txt

#create empty list
selection_list = []

#append to selection list "Name = Project"-
for line in datasetslist:
    selection_list.append("Name = "+"'"+ line+"'")

#from selection list create a sql statement/string '"Name = Project or Name =...."
selection_string = ' or '.join([str(elem) for elem in selection_list])

#use sql to select from a layer
footprint = os.path.join(folder,gdb,"footprint")
footprint_select = arcpy.SelectLayerByAttribute_management(footprint, "NEW_SELECTION", selection_string)

#calculate on selection
arcpy.CalculateField_management(footprint_select,"TITLE",title,"PYTHON")

I assume I can't use the variable as an argument in the SelectLayerbyAttribute_management, however I'm not sure of any other approach. Options?


texts: Is it possible to lock layer order in Layer Browser Panel in QGIS?

texts:
I want to share my projects with my colleauges. The used projects are really huge because of many implemented layers. In order to prevent struggles (wrong overlapping etc.) within the layer order structure in the browser, I am looking for an option to lock layers (or groups) permanently in the browser window.

Is there any way to do this? I could'nt find any suitable tool or script to solve this problem.


texts:
from grass.pygrass.modules.shortcuts import raster as r

gscript.run_command('g.region', raster='C:\\Users\\rit\\Documents\\PhD\\Projects\\Intensity\\input\\2008_cdls\\MI\\cdl_30m_r_mi_2008_albers.tif')

results in the foll. error:

*** CalledModuleError: Module run None ['g.region', 'raster=C:\\Users\\rit\\Documents\\PhD\\Projects\\Intensity\\input\\2008_cdls\\MI\\cdl_30m_r_mi_2008_albers.tif'] ended with error
ERROR: Raster map
Process ended with non-zero return code 1. See errors in the (error) output.
       <C:\Users\rit\Documents\PhD\Projects\Intensity\input\2008_cdls\MI\cdl_30m_r_mi_2008_albers.tif>
       not found

What could be going wrong? the file I am trying to run the command on does exist

-- EDIT: Based on Micha's response:

import grass.script as gscript
import grass.script.setup as gsetup
from grass.pygrass.modules.shortcuts import raster as r

gscript.run_command('r.in.gdal', input="C:\Users\rit\Documents\PhD\Projects\Intensity\input\2008_cdls\MI\\cdl_30m_r_mi_2008_albers.tif", flags='e',output='C:\Users\rit\Documents\PhD\Projects\Intensity\input\2008_cdls\MI\cdl_30m_r_mi_2008_albers')
gscript.run_command('g.region', raster='C:\\Users\\rit\\Documents\\PhD\\Projects\\Intensity\\input\\2008_cdls\\MI\\cdl_30m_r_mi_2008_albers.tif')

I still get this error:

ERROR: Unable to open datasource
it\Documents\PhD\Projects\Intensity\input�8_cdls\MI\cdl_30m_r_mi_2008_albers.tif>
Traceback (most recent call last):
  File "C:/Users/rit/Documents/PhD/Projects/GLAM/Code/Code/Conditions/misc_codes.py", line 88, in <module>
    gscript.run_command('r.in.gdal', input="C:\Users\rit\Documents\PhD\Projects\Intensity\input\2008_cdls\MI\\cdl_30m_r_mi_2008_albers.tif", flags='e',output='C:\Users\rit\Documents\PhD\Projects\Intensity\input\2008_cdls\MI\cdl_30m_r_mi_2008_albers')
  File "C:\Program Files\GRASS GIS 7.0.4\etc\python\grass\script\core.py", line 395, in run_command
    return handle_errors(returncode, returncode, args, kwargs)
  File "C:\Program Files\GRASS GIS 7.0.4\etc\python\grass\script\core.py", line 313, in handle_errors
    returncode=returncode)
grass.exceptions.CalledModuleError: Module run None ['r.in.gdal', '-e', 'input=C:\\Users\rit\\Documents\\PhD\\Projects\\Intensity\\input\x808_cdls\\MI\\cdl_30m_r_mi_2008_albers.tif', 'output=C:\\Users\rit\\Documents\\PhD\\Projects\\Intensity\\input\x808_cdls\\MI\\cdl_30m_r_mi_2008_albers'] ended with error
Process ended with non-zero return code 1. See errors in the (error) output.


texts: How to export geometry as a table in QGIS


texts: Relationship between scale, ppi, and map resolution

texts:
Wanting to know the screen resolution (pixels per inch) of my device, I displayed my map resolution and scale using this openlayers code:

console.log("resolution: " + map.getResolution());
console.log("scale: " +map.getScale());

which gave the result below:

"resolution: 611.4962261962891" mapTest.html:119
"scale: 1733376.0653978775"

Then I used this formula (from https://msdn.microsoft.com/en-us/library/aa940990.aspx):

  Map scale = 1 : (ScreenRes pixels/inch * 39.37 inches/meter *
  156543.04 meters/pixel * cos(latitude) / (2 ^ zoomlevel))

With my values:

1733376.065 = (ScreenRes ppi * 39.37 inches/meter * 611.496 meters/pixel * cos(51.02)) / (2^8)

The latitude I used is 51.20. But the formula produced 28889.6. I think this is pixels per meter so I divided by 39.37 to get pixels per inch. But that resulted in 733.797 which is an unrealistic pixels per inch. I was expecting value between 72 to 96.

What is going on here? Did I do something wrong or is the formula is wrong? The image is an OpenStreetMap image and the projection used is 900913. I hope some one can help.


texts: Editing and changing the default symbol of North arrow in QGIS 3.10

texts: The North arrow shortcut in the layout environment is very useful and I wonder if I can edit the default symbol and change it, with another symbol, from SVG library?


texts:
I'm using QGIS 3.10.4. I have a vector layer that contains more than 200 cities as multiple points on the layer. The attribute table has one column: city name.
My goal is to keep only 10 certain cities on the layer and remove the other 190 and put these 10 cities on a new layer.
I can manually select the city names and save the selected features into a new layer.
I was wondering if there's an easier way? For example, write a query with the target city names. Because sometimes I may have about 1000 cities, makes manual selection too time-consuming.


texts: How to calculate a field in using updatecursor?

texts:
I'm trying to simply add a string in front of a group of numbers from another field in a new field.  The filed I'm using update cursor in is GIS_ID_PAM and I'm trying to bring over numbers from a column called PointID.  My intent is to get "PAM - " in front of the numbers from the PointID column.  Here's what I have that isn't working.

        urows = arcpy.UpdateCursor(fc)
        for urow in urows:
            urow.DateImport = str(timestr)
            urow.GIS_ID_PAM = "PAM - " + [PointID]
            urows.updateRow(urow)


texts:
I'm trying to convert a matplotlib contourf to geojson, but i'm stuck.
I created contourf from a csv.

Here's my code :

data = np.genfromtxt('centroid.csv', delimiter=',')
x = data[..., 0]
y = data[..., 1]
z = data[..., 2]
xi = np.linspace(x.min(), x.max(), 200)
yi = np.linspace(y.min(), y.max(), 200)
zi = griddata(x, y, z, xi, yi)

nb_class = 3

collec_poly = plt.contourf(xi, yi, zi, nb_class, vmax=abs(zi).max(), vmin=-abs(zi).max())

geojson = geojsoncontour.contour_to_geojson(contour=collec_poly, contour_levels=nb_class, geojson_filepath='out.geojson', ndigits=3, unit='m')

it returns me this error :

assert len(contour_levels) == len(collections)
TypeError: object of type 'int' has no len()

I can't really understand all needed arguments for contour_to_geojson and I didn't find a good documentation.

Somebody can help me ?


texts: Scaling x y coordinates in qgis print composer while using atlas generator

texts:
So I am using QGIS 2.18.0 and have made a map template that uses the atlas generator in the print composer. It uses the grids settings for the map to add UTM north and Eastings along the sides of the map.

Is it possible to set the grids to adjust based on the scale of the map?

For example, the x and y interval for the grid would decrease automatically when the atlas generator is producing a map at 1:10 000 and then increase when the producing a map at 1:100 000? Similar to how the visibility of different layers can set based on scale.

Right now it appears that the grid can only be set to one interval.


texts:
For an examination of a local power grid, I need to mark the transformation stations (substation) by simply create a point for each station.

What I want to do is that QGIS creates a point at the end of a line automatically.

How can I get there?


texts:
I have two map frames on my composer. One is a regional, large scale map, and the other is a local, small scale map. I would like the scale of the regional map to automatically update based on the scale of the local map.

For instance, my local map is map 0. I would like the expression in my scale field for my regional map (Map 1) to be 'Map 0 scale + 35000'. Is there a way to do this in the composer?


texts: Python script tool fails when processing rainfall data from Google Earth Engine

texts:
I am using the Earth Engine Python API in a script tool that I'm running from ArcMap. Essentially, the tool will allow the user to draw an area of interest, and from this area of interest, the total rainfall will be calculated using CHIRPS data from Google Earth Engine.

I am able to get this to work fine when I run each line individually from the IDLE development environment: rainfall is calculated within my area of interest and a bar graph is returned. However, when I run the script from the script tool in ArcMap, I get the following error:

  an error occurred: ImageCollection.getRegion: Too many values: 3452924 points x 1 bands x 11 images > 1048576.

I tried drawing a very small area of interest (ten square meters), but I still I got this error. The problem is occurring when I use the getRegion command to bound the image collection by just the area of interest. But again, this didn't throw an error when I was using IDLE.

Below is the relevant snippet of my script (note: the getRegion command is in the third-to-last line of the code):

def earthEngine(aoi):
    ee.Initialize()
    arcpy.AddMessage("In GEE function")
    #Set the bounding box coordinates of the aoi
    desc = arcpy.Describe(aoi)
    xmin = desc.extent.XMin
    xmax = desc.extent.XMax
    ymin = desc.extent.YMin
    ymax = desc.extent.YMax
    boundingBox = ee.Geometry.Polygon([[[xmin, ymax], [xmin, ymin], [xmax, ymin], [xmax, ymax]]])

    #CHIRPS image collection
    chirps = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY')

    #Set the time period over which images will be captured.
    #Only rainfall from the past year will be reported.
    now = datetime.datetime.now()
    endYr = now.year
    endMo = now.month
    endDy = now.day
    stYr = endYr - 1
    stMo = endMo
    stDy = endDy
    #Create variables for years and months of interest
    years = range(stYr, endYr + 1)
    endMoRange = endMo + 13
    months = range(endMo, endMoRange)
    #Set date in EE format
    startdate = ee.Date.fromYMD(stYr, stMo, 1)
    enddate = ee.Date.fromYMD(endYr, endMo, 1)
    #Filter CHIRPS by date
    Pchirps = chirps.filterDate(startdate, enddate).sort('system:time_start').select('precipitation')
    #Function to calculate monthly precipitation total
    def calcMonthlySum(imageCollection):
        mylist = ee.List([])
        for m in months:
            if m <= 12:
                mo = m
                yr = years[0]
                w = imageCollection.filter(ee.Filter.calendarRange(yr, yr, 'year')).filter(ee.Filter.calendarRange(mo, mo, 'month')).sum()
                mylist = mylist.add(w.set('system:time_start', str(mo) + '/' + str(yr)))
            else:
                mo = m - 12
                yr = years[1]
                w = imageCollection.filter(ee.Filter.calendarRange(yr, yr, 'year')).filter(ee.Filter.calendarRange(mo, mo, 'month')).sum()
                mylist = mylist.add(w.set('system:time_start', str(mo) + '/' + str(yr)))
        return ee.ImageCollection.fromImages(mylist)

    #Call function to calculate monthly rainfall
    monthlyChirps = ee.ImageCollection(calcMonthlySum(Pchirps))
    #Limit output to the area of interest
    monthlyChirps2 = monthlyChirps.getRegion(boundingBox,24000,"epsg:4326").getInfo()

    #Create pandas dataframe to plot
    df = pd.DataFrame(monthlyChirps2, columns = monthlyChirps2[0])

aoi, which is passed to the earthEngine function, is a polygon that the user draws.

It seems that the problem is that the data being passed to getRegion is too large. I'm curious if there's a work-around to this, and/or if there's an explanation as to why this will work in an integrated development environment, but not when it's run from ArcMap.


texts: Using GPS tracking with QGIS and Panasonic Toughpad

texts: I'm using QGIS with a Panasonic Toughpad FZ-M1 with inbuilt dedicated GPS. I can connect okay to the GPS using View --> Panels --> GPS information. Satellites and signal strength show up, but no GPS information shows up (i.e. no positional data or any other attributes) in the GPS information window, and no cursor shows on the map. The GPS reports positional data okay when connected to GPS Viewer tool, so the problem appears to be specific to QGIS and the Panasonic Toughpad. I have tried the solutions listed here (Using the Live tracking with internal GPS?) but they do not work for me. Does anyone have any other suggestions to fix this problem?


texts: How to use Canvas element with OpenLayers?

texts:
I would like to use element canvas of HTML 5 in OpenLayers to render maps.

I am searching some code or reference in which I can check it, can you help me in this?


texts: Installing LARCH-EcoPlanIt toolbar in ArcMap?

texts:
In my study I´m looking for ecological corridors and I would like to use the LARCH model for it. I have downloaded it from this side http://www.wageningenur.nl/nl/Expertises-Dienstverlening/Onderzoeksinstituten/Alterra/Faciliteiten-Producten/Software-en-modellen/LARCH.htm . Now I would like to add the tool to my ArcMap (version 10.2) document. The instruction manual already says that the tool works best in version 9.2 and 9.3 but I want to give it a try.
To insert it into ArcMap the instruction manual says that I have to unzip it and then to install it by clicking on setup.exe and then to dubbelclick on clusdist.exe. But if I look at the unzipped file, there are no .exe files at all in it.


texts: Programmatically changing layers in QGIS Print Composer?

texts:
I'm trying to programmatically generate new layers from a database, and then print them using a template. Using the QGIS 2.14.10x64 Python Console in Windows 7x64. I've been following on to code from this answer and this one. Set up code:

from qgis.PyQt.QtXml import QDomDocument
from qgis.utils import iface

def load_print_composer(template):
    '''Load a print composer template from provided filename argument

    Args:
        template: readable .qpt template filename
        console: boolean if method is used in QGIS console

    Returns:
        myComposition: a QgsComposerView loaded from the provided template
        --mapSettings: a QgsMapSettings object associated with myComposition'''
    # Load template from filename
    myDocument = QDomDocument()
    with open(template, 'r') as templateFile:
        myTemplateContent = templateFile.read()
        myDocument.setContent(myTemplateContent)
    composerView = None

    composerView = iface.createNewComposer()
    composerView.composition().loadFromTemplate(myDocument)
    myComposition = composerView.composition()
    mapSettings = myComposition.mapSettings()
    return {'QgsComposition': myComposition,
            'QgsMapSettings': mapSettings,
            'QgsComposerView': composerView}

template = 'somefilename'
printcomposer = load_print_composer(template)

Things that don't work

#1 changing QgsMapSettings:

layerlist = printcomposer['QgsMapSettings'].layers()
layer0 = layerlist.pop(0)
printcomposer['QgsMapSettings'].setLayers(layerlist)

#2 Making sure mapCanvases are linked

From this answer

map0 = printcomposer['QgsComposition'].getComposerMapById(0)
map0.renderModeUpdateCachedImage()
map0.setMapCanvas(iface.mapCanvas())
map0.layerSet()
map0.setLayerSet(layerlist)


texts: How do I transform CSV data into Shapefile in ArcGIS Pro ?


texts: Is transformation of spatial data lossy?

texts:
I have a shapefile that contains data of a country from the following Geographic Coordinate System ( GCS_North_American_1927 aka NAD27 )

contents from the shape file .prj:
GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

I have the following questions regarding this data.

1.) When i load this shapefile into SQL server in NAD27 Geographic coordinate system and also convert it to WGS 84 and load it into sql server.  When i render both on google maps they look as if they overlap perfectly.  Is this because google doesn't have the precision to go out to the proper decimal places to render the spatial data?

NOTE: Google maps assumes you will give it WGS 84 spatial data to work with.

2.) Are Transformations lossy? for example if i perform the following transformation

NAD27 -> WGS 84
and then go back
WGS -> NAD27

do i get the same exact spatial information back or is some precision lost?

3.) If i transform my information from NAD27 to WGS84 to be displayed on google maps will my coordinate alway show the right position of the spatial data relative to WGS84 or will it be slightly off?


texts: Calculating percentage of area from polygons inside one polygon in ArcMap?

texts:
I need to calculate percentage of different areas contained in one Cadastral parcel.
Every parcel (polygon) is overlapped with one or more polygons. I need to find percentage of every overlapping polygon inside parcel. I have used intersect command to get intersected polygons (attribute table on image below).

How can I group values by one attribute so I can later calculate percentage of areas in another field?

How can I group values by one attribute (BROJ_KC) to get something like this (green part of table)?:


texts: How to clip a big shapefile from another shapefile's extent

texts:
I am working on mangroves in Sundarbans (Bangladesh / India), using Landsat pictures.
I need to equalize them in reflectance, before classification, for the evolution of the biomass in this region.
I have a mangroves' shapefile over this area, with +150 000 features (polygons).
So far, the segmentation part and the classification part work well. However, to be able to equalize the pictures, I need to focus on the overlapping parts between this pictures.
I succeeded in extracting the overlap between two pictures, and then having the shapefile corresponding to the boundaries of the overlap.

First, what I need to do is clipping the big mangroves' shapefile to the extent of an overlap between 2 images, so I can then compute the overlapping raster with only the mangroves parts to equalize these two pictures.
I tried this following command :

ogr2ogr -f ESRI Shapefile -clipsrc extent.shp mangrove_clip.shp mangrove.shp

I am coding in python, but even on a direct command line request, this doesn't seem to work. Actually, I don't have any message to screen here because nothing is screened on the terminal, as if the command ran indefinitely...
I know that on QGIS, there is the clip tool under "Vector > Geoprocessing Tools". This works perfectly, and it is really fast (less than 30s).
The thing is that I need to script this because I have over 10 images on this area and I need to work on every overlap between every images, so doing it manually through QGIS is not feasible. Besides, I may have to re-use that for other parts of the world.

That's a pity that for the vector part, QGIS does not screen the corresponding command line like it does for the raster (filling the gdalwarp with its arguments for example, when you want to clip a raster).
Does anybody know what's wrong with my ogr2ogr command ? i know that's a hard question considering I have no returned error message...
What's the command QGIS is using ?
I tried hard to find an answer but every topic ended with this ogr2ogr command as a solution.

And the command actually works for another shapefile...

Even if by letting the command run for 2 days, it worked, I don't think it's normal for it
to take so long when QGIS does it in 30s.

Concerning the mangroves' shapefile (which looks to be the problem, but once again, it works with QGIS), it is a polygon shapefile, and my extent shapefile is also a polygon shapefile.


texts:
I am working on fusion of two different rasters. One raster is a multi-spectral band having spatial reference GCS-WGS-1984 and other raster is SAR band having spatial reference WGS-1984-polyconic. Prior to fusion, I have registered multi-spectral band with SAR band. So spatial reference of multi-spectral data has been changed to WGS-1984-polyconic. After fusion, I want to compare fused raster with original multi-spectral raster. I tried for re-projecting fused raster to GCS-WGS-1984 in QGIS and ArcMap. But during this process, size of my raster is changed from size of original multi-spectral raster.

How do I overcome this problem in QGIS?


texts: How to warp 3D model to take into account horizon curvature?

texts:
I am using USGS Astrogeology data for the Moon's south pole. Here's the same dataset at half resolution:
"LRO_LOLA_DEM_SPolar875_10m.tif (tif) 278 MB" on this page.

I already have a process for converting the GeoTIFF into a 3D model, but because the GeoTIFF is a map of planet radii, it doesn't include any horizon curvature. This results in the outer extents of it being more than a kilometer above the actual terrain. How should I correct for this?

QGIS GeoTIFF info:


texts: gdal_translate: error 1: Could not find an appropriate zoom level

texts:
I'm trying to put a jpeg file (a map) into Sqlite mbtiles.

The command:

gdal_translate -gcp 0 0 6.756394 45.192361 -gcp 0 1080 6.756394 45.127057 -gcp 1820 0 6.912606 45.192361 -gcp 1820 1080 6.912606 45.127057  etache.jpg etache.mbtiles

returns the following error:

Input file size is 1920, 1080
ERROR 1: Could not find an appropriate zoom level

Any clue?


texts: Authentication using PG 10 scram-sha-256 password encryption

texts:
My PostgreSQL 10 database is configured to encrypt passwords using the scram-sha-256 algorithm. I can't establish a connection from QGIS 2.18, as it fails with the following error:

  authentication method 10 not supported

What can I do to connect to the database using the scram-sha-256 password encryption?


texts:
I am using OpenLayers 2 for making clustering on vector point data for which I am a beginner. I have installed GeoServer locally via opengeo suite and using 'world_cities' dataset as an example. I have been able to implement  simple cluster strategy. Now I want to display count numbers on each cluster. For example, when I start to zoom in , the big cluster should display how many element points  are inside. Then when I zoom in further, the big cluster will break down and at the same time show  it shoud display the reduced number of points and until it becomes 1.

Here is my code.

At which position shall I place another method/rule to get numbering staff done ?

var style = new OpenLayers.Style({
                pointRadius: "${radius}",
                fillColor: "#FF3300",
                fillOpacity: 0.8,
                strokeColor: "#66FF00",
                strokeWidth: "${width}",
                strokeOpacity: 0.8
            }, {
                context: {
                    width: function(feature) {
                        return (feature.cluster) ? 2 : 1;
                    },
                    radius: function(feature) {
                        var pix = 2;
                        if(feature.cluster) {
                            pix = Math.min(feature.attributes.count, 7) + 2;
                        }
                        return pix;
                    }
                }
            });

  var wfs_layer2 = new OpenLayers.Layer.Vector("world_cities", {
  strategies: [new OpenLayers.Strategy.BBOX()],
  protocol: new OpenLayers.Protocol.WFS({
  version: "1.1.0",
  url: "http://localhost:8081/geoserver/wfs",
  featureNS :  "http://world.opengeo.org",
  featureType: "cities",
 }),

   strategies: [
            new OpenLayers.Strategy.Fixed(),
            new OpenLayers.Strategy.Cluster({distance: 30})
            ],
            styleMap: new OpenLayers.StyleMap({
                    "default": style,
                    "select": {
                        fillColor: "#8aeeef",
                        strokeColor: "#32a8a9"
                    }
                })

});


texts:
This may be simple question but I can't get started. This tutorial from pen state university has xml request

<wfs:GetFeature service="WFS"  version="1.0.0"
  outputFormat="GML2"  – we would like gml2 as the return format
  xmlns:topp="http://www.openplans.org/topp"
  xmlns:wfs="http://www.opengis.net/wfs"
  xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/wfs
  http://schemas.opengis.net/wfs/1.0.0/WFS-basic.xsd">
  <wfs:Query typeName="topp:bc_roads">  we want to query
                bc_roads features
    <ogc:Filter>
      <ogc:FeatureId fid="bc_roads.1"/> only return features
                                        with an ID of bc_roads.1
    </ogc:Filter>
    </wfs:Query>
</wfs:GetFeature>

My question here is how do I send this request to geoserver using geoext or javascript?


texts: Could not load file or assembly : ESRI.ArcGIS.ADF.Local.dll

texts:
i worked in an application vb.net , arcengine : i copy the project source and i paste it in other Machine , Ooops i can't compile it , Error is :

Error 1 Could not load file or assembly 'file :/ / / C :/ Program Files (x86) / ArcGIS/DeveloperKit10.0/DotNet/ESRI.ArcGIS.ADF.Local.dll' or one of its dependencies. Attempt to load a program with an incorrect format. Line 142, position 5. C: \ Users \ m\ Documents \ Visual Studio 2010 \ Projects \projetc1\reco.resx 142 5 recM

any suggestion ?


texts: Optimizing osm2pgsql imports for OSM data

texts:
I'm currently building an instance on EC2 on which to import the entire Planet.osm snapshot of the whole Earth's worth of data for some projects we're working on.  I've spun up a large Ubuntu x64 instance and attached plenty of separate storage on an EBS volume for the Postgres database and modified it to house the PGSQL data there.

Now the server is having trouble using osm2pgsql to import the snapshot...  After a couple of attempts with different memory configs and whatnot, the process keeps outputting "Killed" after getting most of the way through; once it was killed while "going over pending ways" and the next time, after slightly adjusting the slim cache, it reached "processing ways" before crashing out.  From what I've read, this is generally due to memory issues.

Here's my latest attempt to run the import:

osm2pgsql -v -U osm -s -C 4096 -S default.style -d osm /data/osm/planet-latest.osm.bz2

And here are the specs for a Large instance on EC2:

  Large Instance 7.5 GB of memory, 4 EC2 Compute Units (2 virtual cores with 2 EC2 Compute Units each), 850 GB of local instance storage, 64-bit platform

My question is -- are there some good benchmark resources to determine the tuning requirements for osm2pgsql and Postgres?  Speed of import isn't even that important to me, I'd just like to be able to make sure the process completes safely, even if it takes 4 or 5 days...  I've read through Frederick Ramm's "Optimising the rendering chain" (PDF) document from last year's SOTM, but are there other good opinions / resources?


texts: Does anyone know any sources for free European city boundary data?


texts: Arcmap - add an image to a legend that is not part of layer symbology

texts:
I am making a simple map, showing store locations on a background. Each store brand is displayed with an image (a coloured square with a letter e.g. Aldi might be a blue square with the letter "A" inside).

In the legend we want to display: Brand logo image - Symbology image - Name.
Something like this:

How do I add the brand logo image to the legend?

I don't want to convert the legend to graphics, because I'm exporting a number of maps displaying different layers. I'm thinking of doing the following:

Create a combined image of brand logo and symbol (there's only half a dozen of these so it's manageable)
Duplicate the layer, but "move" the coordinates to display outside the visible map extent
Use the combined images for the new layer symbology
add the layer to the legend
Remove the actual layer from the legend
(The reason for the last two steps are so that I can still have the legend only showing switched on layers)

I think that would work but it seems inelegant.

Is there a better way of doing this?

E.g. is it possible to add an HTML tag to the legend item to add the logo?


texts: How to use GEOS/C++ to efficiently find all point pairs closer than a threshold?

texts:
Can anyone tell me how to use GEOS/C++ to efficiently find all point pairs in a dataset closer than a threshold distance d?

I suspect this might involve quadtrees or rtrees but not quite sure how to query them.


texts:
I am a complete beginner with ArcGIS and I have been looking for a way to live feed tweets to ArcGIS desktop. A quick google search brings this blogpost up from 2009: http://blogs.esri.com/esri/arcgis/2009/12/15/arcgis-explorer-and-geo-service-integration-twitter-add-in/

However, the link to the add-on is dead. Does this add-on still exist? If so, does anyone know where. If not, does anyone know why?


texts: GeoPandas causing issues with PyQGIS standalone

texts:
I have a code where I would like to be able to use both GeoPandas and PyQGIS. However, it seems that including GeoPandas is causing issues in the output of PyQGIS.

Here is a small illustration using NYC OpenData school districts and census blocks

import geopandas
import os
import sys

# Add to Python Path
sys.path.append('/Applications/QGis.app/Contents/Resources/python/')
sys.path.append('/Applications/QGis.app/Contents/Resources/python/plugins')

# Initialize QGIS
from qgis.core import *
from qgis.gui import *
qgs = QgsApplication([], False)
QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/MacOS", True)
qgs.initQgis()

# Setup Processing
from PyQt4 import QtCore, QtGui
import processing
  # create dummy iface
class DummyInterface(object):
    def __init__(self):
        self.destCrs = None
    def __getattr__(self, *args, **kwargs):
        def dummy(*args, **kwargs):
            return DummyInterface()
        return dummy
    def __iter__(self):
        return self
    def next(self):
        raise StopIteration
    def layers(self):
        # simulate iface.legendInterface().layers()
        return QgsMapLayerRegistry.instance().mapLayers().values()
iface = DummyInterface()
  # initalise processing plugin with dummy iface object
plugin = processing.classFactory(iface)

# Set paths
os.chdir('path/to/example')
districts = 'School Districts/geo_export_cdfc9016-ca65-4c45-8423-7acc709a42a2.shp'
blocks = '2010 Census Blocks/geo_export_ae8a0155-ccdc-4a7b-a42e-6707282035cc.shp'
district32 = 'School Districts/district32.shp'
outfile = 'clipped_blocks.shp'

# Select school district 32 and save as new file
districts = QgsVectorLayer(districts,'districts',"ogr")
query = '"school_dis"=32'
selection = districts.getFeatures(QgsFeatureRequest().setFilterExpression(query))
districts.setSelectedFeatures([k.id() for k in selection])
_writer = QgsVectorFileWriter.writeAsVectorFormat(districts,district32,"utf-8",None,"ESRI Shapefile",True)

# Clip blocks with school district 32
processing.runalg("qgis:clip",blocks,district32,outfile)

qgs.exitQgis()

If I run this code without the import geopandas, then I get this clipping:

which is correct. But if I do it with the import geopandas, then I get this incorrect clipping:

It seems that simply by importing GeoPandas, it messes up the QGIS processing algorithms. This keeps happening in my code. No errors are thrown, and it's sort of random which processes work and which ones don't. I've had issues with "qgis:clip" and "qgis:dissolve", but other ones like "qgis:deleteduplicategeometries" and "grass7:v.clean" seem to work fine.

Is there an issue with the way I've setup the PyQGIS processing or are these two just fundamentally incompatible? Ideally I'd like to be able to use them in the same code as sometimes it's easier to do a task in GeoPandas than PyQGIS and vice versa.

I am running this on a Mac (OS ver. 10.13) with QGIS 2.18.15, Python 2.17.14 and GeoPandas 0.3.0.


texts: Convert EPSG:4326 to EPSG:3857 in OpenLayers and React

texts: Is there a way to convert coordinates from EPSG:4326 to EPSG:3857 in OpenLayers. Or does any one have a method that does this in Javascript?


texts: Merge all shp files in a folder into one with a new field populated with the source filename

texts:
Want to merge all shp files in a folder into one with a new field populated with the source filename. Not getting any errors but not getting any outputs.

Code originally taken from (How to add field with filename when merging shapefiles with ogr2ogr?)

# merge_shps.py
import os
path = "C:\TEMP\SWS_SUDS\Cadastral_Maps"  # path to your folder of .shp files
merge = "SWS_INSPIRE"                         # this will be the name of your merged result
directory = os.listdir(path)
count = 0
for filename in directory:
    print "\n" + filename + "\n"
    if ".SHP" in filename.upper() and not ".XML" in filename.upper():
        # On the first pass, create a clone and add the filename column.
        if count == 0:
             # Make a clone (matt wilkie)..
             cmd = 'ogr2ogr ' + path + '\\' + merge + '.shp ' + path + '\\' + filename + ' -where "FID < 0"'
             print "\n" + cmd + "\n"
             os.system(cmd)
             # Add the field (j03lar50n)..
             cmd = 'ogrinfo ' + path + '\\' + merge + '.shp -sql "ALTER TABLE ' + merge + ' ADD COLUMN filename character(50)"'
             print "\n" + cmd + "\n"
             os.system(cmd)
        # Now populate the data (capooti)..
        print "Merging: " + str(filename)
        # You'll need the filename without the .shp extension for the OGR_SQL..
        filenameNoExt = filename.replace(".shp","")
        cmd = 'ogr2ogr -f "esri shapefile" -update -append "' + \
         path + '\\' + merge + '.shp" "' + \
         path + '\\' + filename +'"' \
         ' -sql "SELECT \'' + filename + '\' AS filename, * FROM ' + filenameNoExt + '"'
        # Uncomment this line to spit the ogr2ogr sentence to the terminal..
        print "\n" + cmd + "\n"
        os.system(cmd)
        count += 1


texts: Assigning elevation to reclassified raster in ArcScene or ArcMap

texts:
How it it possible to make a 3D view of reclassifed raster map in ArcScene or ArcMap?

I have identified a potenial avalanche area using DEM in ArcMap and now I want to have this map in a 3D view in ArcScene but i can't since it is reclassified from 1 to 4 means I don't know the elevation value in the areas but want to present it in 3D view.

I have attached the image of reclassed file.


texts: Howto get matrixset parameters from QGIS WMTS Layer metadata in Python

texts:
I would love go get matrixSet attribute along with bounding box coordinates and maybe URL from QGIS WMTS Layer.
Matrixset should be possible via
parse_qs(layer.source())["matrixSet"]

But I have no clue, wheter/how to get the bounding box coordinates for the matrix set (along with the URL).
The most informative and closest data source would be
layer.dataProvider().htmlMetadata()

Which contains exactly the type of required information. Browsing the C++ code (providers/qgswmsprovider.cpp), did not enlighten me, how the HTML values are obtained.
Thanks for any hint.


texts:
How do I crop the main layer in which the base map and other features are placed on?
I want to work on a select portion of the map as opposed to the entire layout.
I've used the clip function but this is used for features only.


texts: QGIS field calculator giving area in acres in layer / project that is in UTM

texts:
From all the documentation I've been able to find, using $area in the field calculator should return the area in the layer's projection's units.  I have a layer in UTM (the project is in the same UTM CRS, which uses meters for its units) but $area is giving me acres.  I have double checked this in ArcMap, I used the Calculate Geometry feature, selected acres and arrived at the same value.

It is handy because I was planning on converting to acres myself, but I want to make sure I get consistent results down the road.  I do have acres selected as the unit to be used by QGIS' measure tool.  Does that selection effect the output of $area now?


texts: ArcGIS Desktop Professional Certification - opinions now that it has been a around a few years

texts:
When I first heard about the ArcGIS Desktop Professional Certification, I was interested in obtaining this.  I took the practice exam but did not do too well since there were many questions pertaining to aspects of GIS I simply do not use.  I saw that there was a 2-day online course for $1k to bring me up to speed but I did not that find economical.  I see that there is a study guide book for the Associate Desktop but not the Professional Desktop Certification.

I decided to take the wait-and-see approach and now that it is about 3 years later, I would like to know what people think of this certification.  Those of you that have received it, was it worth it?  How did you prepare for it?  Did you take the course or did you have enough knowledge to pass on your own?


texts: ogr2ogr filegdb to postgis change geometry column and SRID

texts:
I have figured out (through the use of this thread) on how to go from a file GDB to posgreSQL using ogr2ogr. Now what I want to do is import all the feauture classes but change the geometry column to geom.

ogr2ogr -f "PostgreSQL" PG:"host=localhost port=5432 dbname=freelance user=postgres" C:\pathto\oakland.gdb -overwrite -progress

when i run this cmd ogr2ogr command it correctly imports the feature classes into postgresql but the geometry column annoyingly gets renamed to wkb_geometry. I believe there is a way to alter a column for an individual layer but how do you do it for all feature classes?

although I know its recommended to ask this in a new question but if I wanted to project/transform to another SRID while importing into PostgreSQL how would I accomplish this


texts: 3D visualisation alternative plugins in QGIS 3.0

texts: As the qgis2threejs plugin is not available in QGIS3, what are the 3D visualisation alternatives to plugin Qgis2threejs in QGIS 3.0?


texts: Finding points which fall 20 miles or beyond city boundary using MapInfo?

texts:
I have a table of geocoded pupils and a layer with a polygon representing my city boundary. I need to find how many pupils are located 20 miles or more outside of the city boundary.

I've been reading and i think this could possible be achieved by creating a buffer of 20 miles around the city boundary polygon and then using the SQL Select to query the data to find the points within the polygon which would allow me to determine which pupils fall outside of it.

I have calculated points within a polygon before using this method but I've never created a buffer. I would like to know how i can create the type of buffer which would make a polygon which is the same shape as my city but 20 miles distance from each of the nodes.

Is this an approach i should be taking or should i be doing this a different way?


texts: How to Separate out the disconnected regions of a raster into independent rasters using ArcGIS?

texts:
Please refer to the attached image. The overlay on the map shows a raster file with disconnected circles. I want to split each of those circles into independent rasters in one shot.

Is that even possible? If yes how? If not whats the fastest way to do it?


texts: Finding highest vertex of 3d line with arcgis 10.5 and convert to point feature class?

texts:
I have a series of 3D perpendicular lines running along an embankment. I need to determine the highest point along these perpendicular line. Each perpendicular line has a unique ID so that I can theoretically determine the highest point along each 3d line and join them to create a ridge.

My first problem is how to determine and extract the highest vertex of each 3D line and convert to point. I have 3d Analyst as well as the data reviewer. The Data reviewer can check for z-values but cannot determine the highest.


texts: Difference between a GNSS receiver and GNSS antenna

texts:
I´m relative new using GNSS data and I would like to understand a little more about the instruments that I´m planning to use in my future research.

I was looking information on the UNAVCO website about Geophysical instruments and they mentioned "GNSS receivers" and "GNSS antennas" as different instruments, and actually they have tables that compare some brands and models of each equipment.  Also, I was looking online but I cannot find a clear definition and functions of each one.

Before today I thought that the receiver was inside the antenna, but I´m not sure now.

What are the differences in function between a GNSS Receiver and a GNSS Antenna?


texts: Overwriting layers in a geodatabase maintaining its position in a feature class list

texts:
I'm working with ArcGIS 10. I work on ecological habitat mapping.

I have some features classes inside a geodatabase which represent years of data on a specific area. Each has a reliability field with a value from 0 to 10. What I want to do is getting a map of the most accurate habitat data without overlaps. To do so, I'm looping through the layers of the geodatabase couple by couple. For each couple, I test intersection, and if there is some, depending on the reliability value, I perform an erase task where the layer with the lowest reliability value is erased by the layer with the higher reliability value.

The problem is that each time a layer is erased, I want to replace the previous non erased layer by the new erased layer at the same place in the list of layers. Indeed, I need to reuse the last erased layer if there is a new erase task to perform. But in the code below, I don't succeed at replacing the obsolete layer by the erased one.

I have differents errors :

I have a "same name error" while I authorize overwriting.
If I change the name, the output is at the end of the list and I don't succeed at moving it at the right place in the list.

Here is the code :

import arcpy, os
arcpy.env.overwriteOutput = True
arcpy.env.workspace = arcpy.GetParameterAsText(0)
fcList = arcpy.ListFeatureClasses()
count = 0 #help select the a feature class
cnt = 1 #help select the b feature class
tList = len(fcList) #Count the number of elements in the list
for count in range(tList - 1):
    for cnt in range(cnt,tList):
            desca = arcpy.Describe(fcList[count])
            descb = arcpy.Describe(fcList[cnt])
            namea = desca.baseName # Get the name of the a feature to name the out_poly
            nameb = descb.baseName # Get the name of the b feature to name the out_poly
            arcpy.MakeFeatureLayer_management(os.path.join(arcpy.env.workspace, fcList[count]), "fca")
            arcpy.MakeFeatureLayer_management(os.path.join(arcpy.env.workspace, fcList[cnt]), "fcb")
            arcpy.SelectLayerByAttribute_management("fca","NEW_SELECTION")
            arcpy.SelectLayerByLocation_management("fcb","INTERSECT","fca")

            if arcpy.Describe("fcb").FIDSet: #Check for selection
                    rowa = arcpy.SearchCursor("fca")
                    rowb = arcpy.SearchCursor("fcb")
                    for row in rowa:
                            fiabia = row.fiabilite #get the int value of reliability (all entities of a layer have the same)
                    for rows in rowb:
                            fiabib = rows.fiabilite

                    if fiabia > fiabib:
                            out_poly = nameb
                            fcList[cnt] = arcpy.Erase_analysis("fcb","fca", out_poly) #the new fc is added to the end of the list while it should be at the "count" level of the list.

                            fcList = arcpy.ListFeatureClasses()
                            #arcpy.Delete_management(fcList[cnt])
                            #arcpy.FeatureClassToGeodatabase_conversion("fcb",fcList[cnt])
                            #arcpy.Rename_management(out_poly, fcList[cnt])

                    if fiabia < fiabib:
                            out_poly = namea
                            fcList[count] = arcpy.Erase_analysis("fca","fcb", out_poly)
                            arcpy.Delete_management(fcList[count])
                            fcList = arcpy.ListFeatureClasses()
                            #arcpy.FeatureClassToGeodatabase_conversion("fca",fcList[count])
                            #arcpy.Rename_management(out_poly, fcList[count])

                    else:
                            count = count + 1
                            cnt = cnt + 1

I don't understand why I can't proper update the fcList. If I try with a random list, for example, I can easily change any element of the list :

rlist = [a,b,c,d]
rlist[0] = z
print rlist
[z,b,c,d]

Why is it different with my layers ?


texts:
I´m trying to use some AML files in ArcMap 10, but I have noticed that this option is not available for this version.

Does anyone know or tried to "translate" the aml sequences to another language for example python? or another application like QGIS or SAGA.


texts: OpenLayers 3 - Source with multiple available projections

texts:
Is there anyway to specify a (Tile)WMS source that supports multiple projections?

For example, my WMS server supports requests for both EPSG:4326 and EPSG:3395. So if the projection is supported in server-side, the tiles are requested with the appropriate projection, else OpenLayers tries to reproject one of the supported projections.


texts: Cache geodjango query with bounding box filter

texts:
Im using geodjango (postgis) to retrieve points contained inside a bounding box (made by leafletjs), it works so nice, but now I want to cache this query.

A tipical query is one like this:

/points?xmax=193.0078125&xmin=-200.7421875&ymax=49.83798245308484&ymin=29.53522956294848

This kind of query will never be cached, is there a way to cache this?, I mean, achieve some cache hits, assuming would appear more/less points in the viewport than would strictly correspond. Some ideas...

removing decimal position in the numbers
passing zoomlevel as a param should be better/mandatory?

Thanks for your thoughts.

Solved:

Use https://pypi.python.org/pypi/django-geojson to make a simple view which uses (x,y,z) as params to do filter your queryset.
Use http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Tile_numbers_to_lon..2Flat._5 to make the filter in the view
Use https://github.com/glenrobertson/leaflet-tilelayer-geojson/ for making a layer in you map which makes your tile request to the previous view.


texts: Getting forced close when selecting attachment on QGIS 3.2.1?

texts:
I added a new column, labelled it "pics" then I go to properties, layer form, change widget type of "pics" to attachment and set default path to the folder location of my images. Now when I go to attributes and click the 3 dots(...) it goes to the location of the folder then it forced close the QGIS application.

Why am I getting forced close?


texts:
I am using QGIS version 3.6 and Time Manager version 3.2.3 to create a temporal map of the US's energy usage over time. I have a polygonal shapefile of the US with each row of the .dbf file containing a state's name, its USPS code, corresponding region (e.g. new england, midwest), land area, etc., which means that there are multiple entries within the shapefile that have the same region (Rhode Island is New England, and so is New Hampshire). My other file is a .csv file that contains energy consumption data for the region of New England annually from 2017 to 2050.
Here's the first few rows of the energy data:

And here's where one can download the shapefile: cb_2017_us_state_5m

To start off, I joined these two layers by region (both have a region column). After changing the style to graduated and classifying the intervals, the map only showed states in New England, which suggests I joined the layers properly. However, when I checked the resulting attribute table, I noticed that the time column only displayed the last year represented in the data, and that all other values from the csv file were set to their values from that year. Here's what I mean:

I suspect that the other time values were not saved at all, which was confirmed when I added the joined layer to time manager--the only values visible were the values from 2050, and moving the time manager's slider did nothing.

According to this link, Time Manager does not support joined layers, but I don't know if this applies to my situation because I do not get any errors, just a lot of missing years.

I have tried a few things to get around this problem:

I duplicated the rows for each year in the energy data so that I had the same number of rows for each year as there are states in New England (i.e. I have six rows with identical data from 2017 to match the six states in New England, six for 2018, etc.). This did not solve the problem.
I opened the .dbf file that accompanied the shapefile and duplicated the data so that there would be one set of states for each year in the energy file. The DBF file thus had many more rows than it originally had, but the number of rows that showed up in the attribute table remained the same. Additionally, the error with Time Manager persisted.
I tried both of these at the same time, which also did not work.

I feel like I could work around this problem if I could import a new attribute table to use in place of the one that comes with the shapefile, but I am unsure how to do that, and it's very possible it wouldn't solve the problem anyway.

Is there a way to join the two layers in a way that is compatible with Time Manager?


texts: At the moment I have selected numerous vectors in a layer using the Select By tool, but what I want to do now is go to the attribute table and only find the selected features that would be highlighted to be grouped together for ease of use. Is this possible?


texts: Buttons (Ok, Cancel) disappearing in "Apri Modulo" from Attribute Table in QGIS

texts:
QGIS 2.0 has changed the visualisation of a form opened from the attribute table: instead of opening the form in a separate window, the form is opened in the same one.
This is not a problem, but all the buttons disappear...

I created my own form with Qt4 and all the logic with PyQt and the majority of the constraints are associated to the Ok button: now my logic is no longer useful.

Is it possible to re-enable these buttons in QGIS 2.0.1?

I copy here part of my code, maybe someone can help me.

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from definizioni import *

myDialog = None
denominazione = None
denominazioneCheck = None

def formOpen(dialog,layerid,featureid):
  global myDialog
  myDialog = dialog

  global denominazione, denominazioneCheck
  denominazione = myDialog.findChild(QLineEdit,"denominazione")
  denominazione.textChanged.connect(denominazione_onTextChanged)
  if (denominazione.text() == "NULL"):
    denominazione.setText("")
  else:
    temp = denominazione.text()
    denominazione.setText("null")
    denominazione.setText(temp)

  buttonBox = myDialog.findChild(QDialogButtonBox,"buttonBox")
  # Disconnect the signal that QGIS has wired up for the dialog to the button box.
  buttonBox.accepted.disconnect(myDialog.accept)
  # Wire up our own signals.
  buttonBox.accepted.connect(validate)
  buttonBox.rejected.connect(myDialog.reject)
  self.connect(buttonBox, SIGNAL("accepted()"), self.accept)

def denominazione_onTextChanged(text):
  global denominazioneCheck

  if not validateStringNotNull(text):
    denominazione.setStyleSheet("background-color: rgba(255, 107, 107, 150);")
    denominazioneCheck = -2
  else:
    if not validateAlpha(text):
      denominazione.setStyleSheet("background-color: rgba(255, 107, 107, 150);")
      denominazioneCheck = -1
    else:
      denominazione.setStyleSheet("")
      denominazioneCheck = 1

# Valido la form
def validate():
  if denominazioneCheck == -2:
    messageBoxErrore(myDialog,"La denominazione di un comune non puo' essere vuota")
  else:
    if denominazioneCheck == -1:
      messageBoxErrore(myDialog,"La denominazione di un comune non puo' contenere numeri")
    else:
      # Return the form as accepted to QGIS.
      myDialog.accept()

I need that the field of the form respect some constraints, as you can see.

In case that constraints are not fulfilled I don't want that the record is stored in the database (I don't want accept() signal is called!).


texts: display shapefile from postgres using arcmap?

texts:
I have uploaded shapefile using postgis extension in postgres but when i tried to display the shapefile in arcmap, it wont display the data in vector.

this is the image of my connection.

i followed the steps from this link ; https://www.youtube.com/watch?v=C0W1ZgU8GkQ
as you can see, when the data from postgres is dragged into table of contant it should display the vector of the data. but mine just wont show up, really hope anyone can help if im missing something.


texts: Speed up pgRouting point_to_id using new  operator in PostGIS 2.0.0

texts:
The point_to_id function, which assigns a node id to a point based on existing nodes and tolerance, is using a query as follows:

SELECT
id, the_geom
FROM vertices_tmp WHERE
    ST_DWithin(_p,the_geom, tolerance)
ORDER BY ST_Distance(_p,the_geom) LIMIT 1;

This example is very similar to the one used to explain the new KNN index based queries. I modified the query like this:

SELECT
id, the_geom
FROM vertices_tmp WHERE
    ST_DWithin(_p,the_geom, tolerance)
ORDER BY _p <-> the_geom LIMIT 1;

However, as the documentation for '<->' says, the index kicks in only if one of the two operands is constant, as in ST_MakePoint(10,10).

Is there any trick to get the index to be used also in a function?


texts:
I am trying to stack four bands into one TIFF, this is the code I'm using:

            red_raster = rasterio.open(r"red.tif")
            red = red_raster.read(1, masked=True)

            nir_raster = rasterio.open(r"nir.tif")
            nir = nir_raster.read(1, masked=True)

            green_raster = rasterio.open(r"green.tif")
            green = green_raster.read(1, masked=True)

            blue_raster = rasterio.open(r"blue.tif")
            blue = blue_raster.read(1, masked=True)

            out_img = "stack.tif"

            out_meta = red_raster.meta.copy()
            out_meta.update({"count": 4,
                             "nodata": -10000})

            stack = np.dstack((blue, green, red, nir))
            with rasterio.open(out_img, "w", **out_meta) as dest:
               dest.write(stack, 4)

I am getting the error:

ValueError: Source shape (1, 10980, 10980, 4) is inconsistent with given indexes 1

I tried following the steps here: Using Rasterio or GDAL to stack multiple bands without using subprocess commands

but I don't understand exactly what they are doing and am still getting errors.


texts: Extract/Clip LAS (LiDAR lasd) based on shapefile, using ArcGIS Pro and ArcMap 10.6

texts:
I have a few large LAS files which span a big area. I have created an .lasd dataset in ArcGIS Pro and ArcMap 10.6. I want to use a shapefile to define the extent of the LAS file, and extract/filter out the las points that lie within the shapefile.

I have used the Extract LAS tool (tried both ArcGIS Pro and ArcMap 10.6) to do this, but both gave a 999999 error executing function. Have also tried running ArcGIS Pro/ArcMap as admin, etc, but the problem persists.

Does anyone have any ideas?


texts: Find length and width of regular/irregular polygons

texts:
How can I find the length and width of regular and irregular polygons on QGIS? There's > 1000 so can't measure manually with the line tool.

I have perimeter and area but need the other two measurements too.


texts:
I setup wms service using mapserver, the map file is attached . I have created index page, the code is

<!DOCTYPE html>

    <meta charset="utf-8" />
    <title>SOLAP</title>
    <link rel="stylesheet" type="text/css" href="style.css" />

     <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.6.4/leaflet.css" />

The images are repeating as shown in the photo, I tested this web service from first and works without problems.

thanks

The map file, index.html and image output map are here

www.dapboyaca.gov.co/tools/Mapa.zip


texts: Finding percentage of feature in raster using ArcGIS Desktop?

texts:
I am attempting to find the percentage of polygon features at a 5 meter resolution that are located inside each 25km raster cell. I am looking to assign percentages to each cell based on the percentage of a feature appearing in the cell (10%, 20%...).

Is there a tool that can do this procedure?

The raster being limited to only a single value field has prevented me from being able to run attribute queries.

I am using ArcGIS Desktop.


texts: GIST index not used by query optimizer on geometry field

texts:
I'm trying to add elevation curves on an OpenStreetMap tile server I'm currently building.
I created a table in PostgreSQL to store the shapes of the whole world.
Structure was defined by gdal_contour tool as follows :

                                                              Table "public.contour100"
    Column    |           Type            | Collation | Nullable |                   Default                   | Storage | Stats target | Description
--------------+---------------------------+-----------+----------+---------------------------------------------+---------+--------------+-------------
 ogc_fid      | integer                   |           | not null | nextval('contour100_ogc_fid_seq'::regclass) | plain   |              |
 id           | numeric(8,0)              |           |          |                                             | main    |              |
 height       | numeric(12,3)             |           |          |                                             | main    |              |
 wkb_geometry | geometry(LineString,3857) |           |          |                                             | main    |              |
Indexes:
    "contour100_pkey" PRIMARY KEY, btree (ogc_fid)
    "contour100_wkb_geometry_geom_idx" gist (wkb_geometry)

This table contains 45 341 726 rows.

My CartoCSS YAML project references this table like so :

- id: contour100
    geometry: linestring
    Datasource:
      type: "postgis"
      dbname: "myDatabase"
      key_field: ""
      geometry_field: "wkb_geometry"
      extent: "-20037508,-20037508,20037508,20037508"
      table: (SELECT
            ST_Transform(wkb_geometry, 3857) AS wkb_geometry, height
            FROM contour100
            WHERE wkb_geometry && !bbox!
          ) AS c100

When rendering tiles, Mapnik generates queries that look like this one :

SELECT ST_Transform(wkb_geometry, 3857) AS wkb_geometry, height
FROM contour100
WHERE wkb_geometry && ST_SetSRID('BOX3D(-3.402823466385289e+38 -3.402823466385289e+38,3.402823466385289e+38 3.402823466385289e+38)'::box3d, 3857)

This works and renders my elevation curves nicely, but takes several minutes to process the data.

Using EXPLAIN, PostgreSQL outputs this :

Seq Scan on contour100  (cost=0.00..4827531.68 rows=44987512 width=37)
   Filter: (wkb_geometry && '01030000A011[...]'::geometry)

As you can see, I have a spatial GIST index on my wkb_geometry column, but the query optimizer chooses a sequential scan.
Being unfamiliar with geography indexing, I don't understand the reason why a sequential scan is performed.

I've run VACUUM ANALYSE on this table, and tried CLUSTER on the wkb_geometry column as well, but the GIST index was never used.


texts: Translating by adding coordinate constants in ArcGIS Desktop

texts:
Can you add coordinate constants to ArcGIS shapefiles?

I can do this in QGIS with the Translate tool but can't find a similar tool in ArcGIS Desktop 10.6.


texts: Creating balloon for dynamic KML file using OpenLayers

texts:
I want to create a balloon for a KML file that dynamically moves.
The example that I am working on is sundials.html.
Let's say that the KML is moving.
How we would get the balloon to move as well?


texts: Unauthenticated sensitive Information Disclosure - ArcGIS REST services

texts:
Recently a security audit has been done on our servers and applications which has ArcGIS enterprise installed on it.

Our Architecture setup :

We have ArcGIS installed on the server, ArcGIS exposes its data via its restful api for its clients which is authenticated with OAuth2 token security (Out of the box ESRI Feature)
We have mobile applications for both iOS/Android and web which are built on top of ESRI SDK’s.some of the features in the app consumes data from ArcGIS rest services for their functioning.

The audit reported a potential vulnerability stating that one of the rest service disclose sensitive informations including email and phone numbers even for non authenticated users.

consider web application, we make use of esri’s proxy files to manage access to our resources for them.
the vulnerability is found for the following feature server query through the rest interface,

/webapp/proxy/proxy.ashx?https://www.site.com/ArcGIS/rest/services/PublicPortal/xxxxxx/FeatureServer/0/query?f=json&where=EMAIL%20IS%20NOT%20NULL&returnGeometry=true&spatialRel=esriSpatialRelIntersects&outFields=*

I understand this as a problem with the ArcGIS rest services, but I’m not sure about this, if some one could, please clarify the following points

1)what could be done to mitigate this problem.

Kindly please help.


texts: Why can't Basemap's readshapefile read the world_adm0 shapefile?

texts:
I'm using the world shapefile that can be found at https://trac.openstreetmap.org/browser/applications/editors/merkaartor-branches/merkaartor-0.14-fixes/share/world_shp/world_adm0.shp?rev=17052.

import pysal as ps
import ogr, osr
import shapely.geometry
from mpl_toolkits.basemap import Basemap

polygons = ps.open('world_adm0.shp')

# from http://gis.stackexchange.com/questions/129745/writing-a-shapefile-using-ogr-from-shapely-geometry-no-feature-added-error
# create the shapefile
driver = ogr.GetDriverByName("Esri Shapefile")
ds = driver.CreateDataSource("outputlocation.shp")
layr1 = ds.CreateLayer('',None, ogr.wkbPolygon)

for polygon in polygons:
    sh_poly = shapely.geometry.asShape(polygon)

    # create ogr geometry
    polylyr = ogr.CreateGeometryFromWkt(sh_poly.wkt)

    # create the field
    layr1.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))

    # Create the feature and set values
    defn = layr1.GetLayerDefn()
    feat = ogr.Feature(defn)
    feat.SetField('id', 1)
    feat.SetGeometry(polylyr)
    layr1.CreateFeature(feat)

# close the shapefile
ds.Destroy()

map = Basemap(width = 2000000, height = 1500000, resolution = 'l', projection = 'aea', lat_1 = -23, lat_2 = -32, lon_0 = 25, lat_0 = -29)
map.readshapefile('outputlocation', '', drawbounds=True)

results in:

ValueError                                Traceback (most recent call last)
<ipython-input-3-92ac4f1982b1> in <module>()
     33
     34 map = Basemap(width = 2000000, height = 1500000, resolution = 'l', projection = 'aea', lat_1 = -23, lat_2 = -32, lon_0 = 25, lat_0 = -29)
---> 35 map.readshapefile('outputlocation', '', drawbounds=True)

/apps/python/2.7.6/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in readshapefile(self, shapefile, name, drawbounds, zorder, linewidth, color, antialiased, ax, default_encoding)
   2144         info = (shf.numRecords,shptype,bbox[0:2]+[0.,0.],bbox[2:]+[0.,0.])
   2145         npoly = 0
-> 2146         for shprec in shf.shapeRecords():
   2147             shp = shprec.shape; rec = shprec.record
   2148             npoly = npoly + 1

/apps/python/2.7.6/lib/python2.7/site-packages/mpl_toolkits/basemap/shapefile.pyc in shapeRecords(self)
    541         shapeRecords = []
    542         return [_ShapeRecord(shape=rec[0], record=rec[1]) \
--> 543                                 for rec in zip(self.shapes(), self.records())]
    544
    545 class Writer:

/apps/python/2.7.6/lib/python2.7/site-packages/mpl_toolkits/basemap/shapefile.pyc in records(self)
    513         f.seek(self.__dbfHeaderLength())
    514         for i in range(self.numRecords):
--> 515             r = self.__record()
    516             if r:
    517                 records.append(r)

/apps/python/2.7.6/lib/python2.7/site-packages/mpl_toolkits/basemap/shapefile.pyc in __record(self)
    478                     value = float(value)
    479                 else:
--> 480                     value = int(value)
    481             elif typ == b('D'):
    482                 try:

ValueError: invalid literal for int() with base 10: '**********'

What exactly has gone wrong with writing the polygons?


texts: How can I snap one set of points to another in R

texts:
I would like to snap one set of points (SpatialPoints* objects) to another in R using Euclidean distance. I'm hoping for a function like maptools::snapPointsToLines.

The attributes of the points don't need to be transferred.

Any ideas?


texts: No Option to 'Extract' under 'Raster' tab

texts: So I'm running an earlier version of QGIS (3.4.6, as it is supposed to be more stable I'm told) for a Field Geology course, however after installing it and adding DEM data I am unable to 'Extract>Contours' as the option is not there. Is there a tool I don't have active?


texts: NameError: name 'QgsApplication' is not defined in python standalone script

texts:
I have this standalone script to display shapefile. I followed the qgis python cookbook procedure, but i got this error. "NameError: name 'QgsApplication' is not defined". Here is my whole code. Someone please tell me where am wrong. I verified my setPrefixPath by running it in Qgis python console

from PyQt4.QtGui import *
from PyQt4.QtCore import *
import sys, os

class Mapbox(QMainWindow):
    def __init__(self, shapefile):
        QMainWindow.__init__(self)
        self.setWindowTitle("Map Box Viewer")
        canvas = QgsMapCanvas()
        canvas.useImageToRender(False)
        canvas.setCanvasColor(Qt.white)
        canvas.show()

        layer = QgsVectorLayer(shapefile, "layer1", "ogr")
        if not layer.isValid():
            raise IOError("Invalid shapefile")

        QgsMapLayerRegistry.instance().addMapLayer(layer)
        canvas.setExtent(layer.extent())
        canvas.setLayerSet([QgsMapCanvasLayer(layer)])
        self.setCentralWidget(canvas)

app = QApplication(sys.argv)
QgsApplication.setPrefixPath("C:\OSGeo4W64\apps\qgis", True)
QgsApplication.initQgis()
viewer = Mapbox(r"shapefile.shp")
viewer.show()

app.exec_()
QgsApplication.exitQgis()


texts:
I accept I am being a pedant here, but feel the semantics are important..

If I were to drawn a path representing a 'route' between 3 random points on the earth's surface, then establish the nearest elevation points, say every 20m along that path, would you consider that set of triples to be the:

The 'Route Elevation'
The 'Route Contour'
Something else

The definition seems to be used interchangeably, and I'm not sure which one is (more) correct.

The definition of 'Contour' - as I understand it, has a specific meaning of equivalence ie. the line is of constant value along its length, which a path containing various elevations wouldn't be.

Yet the definition of 'Elevation' is a height above a specific datum point, rather than a collection of triples describing it.

The reason this data is an array of triples not just an array of elevation values is the route can be explicitly declared as a pair of lat/longs, however my elevation data may not fall exactly along that path, and so will be an array of triples separate from the route data.


texts: Error while importing CSV File using PyQGIS

texts:
Currently I am trying to import a CSV file using the python extension in QGIS and I have run into a problem that I can not seem to solve.

I am trying to do a basic import of a CSV file from a python code I created (See below) and every time the file loads the values in my attribute table become re-arranged (see image).

When I load the file from the "Create a Layer from a Delimited Text File" option it loads without a problem.

Here is a sample of my code:

uri = "file:///D:/MLB Stadiums/Test_R2/Player_Files/Espinosa_2014.csv?
       delimiter=%ss&xField=%s&yField=%s" % (",", "field_3", "field_4")
   vlayer = QgsVectorLayer(uri, "Espinosa", "delimitedtext")
   caps =vlayer.dataProvider().capabilities()

Here is a comparison of what it should look like after import and what it does look like.

Does anyone know why it is rearranging the values?


texts: Results disappear after SAGA toolbox operation completes (QGIS)

texts:
I'm running into this problem on QGIS 3.18 and 3.16.
When I run the SAGA multiple regression analysis-tool, there seems to be no issues in completing the task. From the log-window the following appears:
C:\OSGeo4W64\bin>exit
Execution completed in 47.98 seconds
Results:
{'INFO_COEFF': 'C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/f6fe324738e54a47a620ebe3fcf49a1a/INFO_COEFF.shp',
'INFO_MODEL': 'C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/06f08311e82e4a09bc3f1dcbe65cbd83/INFO_MODEL.shp',
'INFO_STEPS': 'C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/82fbd5f3529e43989f8f990deebe05e9/INFO_STEPS.shp',
'REGRESSION': 'C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/b70cf6e190be41e2b1b3b19c96a97e74/REGRESSION.sdat',
'RESIDUALS': 'C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/c1f03da107c04d608de0c84161f6a893/RESIDUALS.shp'}

However, immediately thereafter, I get the following errors:
Loading resulting layers
The following layers were not correctly generated.
• C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/06f08311e82e4a09bc3f1dcbe65cbd83/INFO_MODEL.shp
• C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/72bc589caf8d41638e4d0cd3a2d7c36d/REGRESSION.sdat
• C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/82fbd5f3529e43989f8f990deebe05e9/INFO_STEPS.shp
• C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/b70cf6e190be41e2b1b3b19c96a97e74/REGRESSION.sdat
• C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/c1f03da107c04d608de0c84161f6a893/RESIDUALS.shp
• C:/Users/Gebruiker/AppData/Local/Temp/processing_sASWWc/f6fe324738e54a47a620ebe3fcf49a1a/INFO_COEFF.shp
You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

It is as if the files are deleted instantaneously and can therefor not be loaded. Has anyone ran into this issue before and found a way to work around it?
Things I have tried

Save to a (non-temporary) file at a different location. Result: there are no files there, though they appear to have been there temporarily as I can see (but not open) them from the 'recently opened' section in my file explorer.


texts: Does ArcGIS Server allow the output data to be stored in a database?

texts:
There is a chapter in my bachelor thesis on database output storage support within various WPS implementations. In other words, I want to find out if the output data (data derived from processes executed on a server) can be stored in a database (such as PostGIS).

For example, the ZOO-Project, one of the implementations of the WPS standard, includes such option and explains it on their website:

http://www.zoo-project.org/docs/install/installation.html#zoo-install-db-backend

Now, I need to cover geoprocessing services executed on ArcGIS Server in this chapter as well, but I cannot find the information I need on the Esri website.

So, if I am a consumer of a process executed on an ArcGIS Server, what options do I have as for storing the output data?

Is there an option to store the data within a (spatial) database?


texts: How to determine if a raster layer covers a polygon layer?

texts:
I have a raster and a polygon layer. In my application I always want the raster layer to cover all the polygon layer. How can that be determined? Is it possible to find out from boundary extent of these layers?

In this example, I want to give a warning to the user that the raster layer does not cover the polygon layer.

I am using .Net technology in ArcMap 10.0.


texts: Multiple features within a single feature layer with different geometries

texts:
I am working with ArcObjects SDK 10 for .NET. There is already a geospatial database at my disposal. So I am trying to map its spatial context to ArcObjects context.

We have objects like state, country, city, etc. Each object has one or more geometry definitions that are not necessarily the same type. For example, a state object for California contains the polygon for state shape and additionally a point for the state's center point (centroid).

After reading this article, what I understand is that a FeatureLayer is a UI concept to display layers of Features on ArcMap and FeatureClasses are a way to group same geometry type Features. In addition, there is a one-to-one mapping between a FeatureLayer and a FeatureClass.

So I have two questions:

There is no way a FeatureLayer could contain Features with different types, right?
It is impossible for me to create one layer per my object on ArcMap, because my object has multiple geometry types. Is that right?

Then, how could I come up with a layer definition that lets me select and operate on same object at the same time (such as select cities, select countries)?


texts: ArcGIS 10.2.2 make a specific query with IF

texts:
I'm using ArcGIS 10.2.2 with a geodatabase stored in a file, I'm not using ArcGIS server. I want to do a specific query in it. I want to do a query which look like  :

SELECT *, IF(column1==column2,'1','0') FROM mytable

Is it possible ? Where can I enter a query like that ?


texts:
I am wondering how it really works, I mean "order by" statement when it comes to geometry/geography columns in postgis/postgresql. Besides that really need a way gather rows from a table with a point geometry column with "left" direction. Suppose we have the coastline from an island on a table where each row represents a point on the island coastline, is it possible to select all the points with a query like that?

select * from island_table order by ?


texts: How to prevent "Type error: unhashable type" while writing properties with Fiona?

texts:
Using shapely's dissolve method I created new geometries from two existing shapefiles. Now I would like to save those new geometries and their properties with Fiona.

I have tried to write the properties in Dict as well as OrderdDict format, but still keep receiving Type error: unhashable type: 'Dict' or Type error: unhashable type: 'collections.OrderedDict'.

I consulted several sources, including:

Using Fiona to write a new shapefile from scratch
Fiona - Preffered method for defining a schema
Fiona Manual: https://fiona.readthedocs.io/en/latest/manual.html

1 Data structure

for poly_element in difference_geometry:
    diff_feature = {
        "geometry" : {
            'type': 'Polygon',
            'coordinates': poly_element,
            },

1a - OrderedDict

        # ordered dictionary
        "properties" : OrderedDict([
            ("area_code", difference_area_code),
            ("area_name", difference_area_name),
            ("area_type", difference_area_type),
            ]),
        }

    shapely_geometries.append(diff_feature)

1b - Dict

        # dictionary
        "properties" : {
            "area_code" : difference_area_code,
            "area_name" : difference_area_name,
            "area_type" : difference_area_type,
            },

    shapely_geometries.append(diff_feature)

2 Output parameters

stat_area_crs = from_epsg(28992)
output_driver = "ESRI Shapefile"

2a - OrderedDict

stat_area_schema = {
    'geometry' : 'Polygon',
    'properties' : OrderedDict([
        ('area_code', 'str'),
        ('area_name', 'str'),
        ('area_type', 'str')
    ])
}

2b - Dict

stat_area_schema = {
    'geometry' : 'Polygon',
    'properties' : {
        'area_code' : 'str',
        'area_name' : 'int',
        'area_type' : 'str'
    }
}

3 Write output to shapefile

with fiona.open(
        'new_stat_areas.shp',
        'w',
        driver = output_driver,
        crs = stat_area_crs,
        schema = stat_area_schema) as c:
    for feature in shapely_geometries:
        c.write({
            'geometry': mapping(feature['geometry']['coordinates']),
            'properties': {feature['properties']}
        })

4 Error message

4a - OrderedDict:

Traceback (most recent call last):
  File "dataIntersections.py", line 191, in <module>
    'properties': {feature['properties']}
TypeError: unhashable type: 'collections.OrderedDict

4b - Dict:

Traceback (most recent call last):
  File "dataIntersections.py", line 190, in <module>
    'properties': {feature['properties']}
TypeError: unhashable type: 'dict'

5 Extra info

When replacing the c.write part for this print statement:

pprint.pprint({'geometry': mapping(feature['geometry']['coordinates']),
      'properties': feature['properties']})

The script returns this (cropped) information:

5a - OrderedDict

{'geometry': {'coordinates': (((144200.28399658203, 460830.77520751953),
                               # (...)
                               (150070.26199916803, 454926.5999999999),
                               (150039.87305716722, 454945.22763747524))),
              'type': 'Polygon'},
 'properties': OrderedDict([('area_code', 'GM0355'),
                            ('area_name', 'Zeist'),
                            ('area_type', 'Rural')])}

5b - Dict

{'geometry': {'coordinates': (((153474.06604003906, 469029.67779541016),
                              # (...)
                               (153473.87103271484, 469029.7180175781),
                               (153474.06604003906, 469029.67779541016)),),
              'type': 'Polygon'},
 'properties': {'area_code': 'GM0307',
                'area_name': 'Amersfoort',
                'area_type': 'Rural'}}

I would like to know what goes wrong and how to solve it!


texts: Converting a text (lat/lon) to a point or geography

texts:
In my table I have a column filled with text type lat/lon entries. I want to convert these entries to a point? Or something I can use with PostGIS functions

I am aware of SELECT ST_PointFromText('POINT(-71.064544 42.28787)', 4326);, but it does not work when I try to do SELECT ST_PointFromText('POINT('locdata.loc_geo')', 4326) from locdata

That just gets me this error: ERROR:  syntax error at or near "locdata"
LINE 1: SELECT ST_PointFromText('POINT('locdata.loc_geo')',.

Extra info:
my locdata.loc_geo entries are text entries structured like this: 34.032601, 78.551298 (the order doesnt matter, i dont need an exact place, just get the math working)


texts: Load an image into QGIS by clicking on a vector footprint

texts:
I have a shapefile with image footprints. Each footprint/shape contains links (both local and web) to the georeferenced image it represents.
Is there any way to load an image into QGIS when user clicks on a corresponding shape?


texts: Issue in calculating NDVI using Rasterio Python

texts:
I am trying to calculate NDVI using two clipped raster images of Landsat 7 (NIR & Red Bands clipped using mask file) using the following code:
import rasterio as rio
import numpy as np
import matplotlib.pyplot as plt

with rio.open(r'D:\clip_test_b3.tif') as src:
    red = src.read(1)   # (Rows, Columns) = (2731, 3660)

with rio.open(r'D:\clip_test_b4.tif') as src:
    nir = src.read(1)   # (Rows, Columns) = (2730, 3635)

np.seterr(divide = 'ignore', invalid = 'ignore')

ndvi = (nir.astype(float) - red.astype(float))/(nir + red)

plt.imshow(ndvi)

In the above code both the bands (Red & NIR) are of different shapes (different rows and columns).
After running the above code I am getting message "ValueError: operands could not be broadcast together with shapes (2730,3635) (2731,3660) ".
But when same NDVI calculation I am trying to do in ArcMap (using Raster Calculator), then NDVI is getting calculated.
Can someone please help me out in solving out this error.