texts:
I am new to Python and scripting in general.

My overall goal is to create a script that I can cycle through all the raster maps in a mapset and then save the univariate stats to a .csv file for each raster map.  I am working on the first part: cycling through all the maps and printing the stats to the terminal window. Later, I will try and save the data to a .csv file.  I have seven maps in my practice mapset and I get seven errors, but my problem is that I do not know what variable to assign to the parameter MAP (hopefully, I stated that correctly).

Here is what I have thus far:

import sys
import os
import grass.script as grass
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb="/Users/kc/GrassData"
location="Data_2013"
mapset="methods_check"
gsetup.init(gisbase, gisdb, location, mapset)
grass.run_command("g.list", _type="rast")
mymaps = grass.parse_command("g.list", _type="rast")
print len(mymaps)
for items in mymaps:
    grass.run_command('r.univar', map='items', separator=",")

What should I put for map = 'items' so that it cycles through all the raster maps of a mapset?

I am using GRASS 7.0, but the same thing happens in GRASS 6.5.  In Grass 6.5, the parameter separator must be changed to fs. Both versions label the input file "map".


texts: Using non-unique names for attribute indexes of file/ArcSDE geodatabases?

texts:
I am using ArcGIS 10.2.2 for Desktop.

Within a file geodatabase does an attribute index name need to be unique?

I thought that it did, but I have not been able to find any documentation to say whether it does or not, for file or any type of geodatabase.

I have tested by:

copying a test feature class from one file geodatabase into another as TestFC
creating an attribute index on a field in TestFC
renaming TestFC to be called TestFC_OLD
copying the same original test feature class in again as TestFC
creating an attribute index on a field in the new TestFC using the same attribute index name
Looking at the Indexes tab of the Properties for both TestFC_OLD and TestFC to see that both have an attribute index of the same name

ArcGIS for Desktop seems to have no objection to the above workflow, so my only concern would be if queries to one feature class were able to use an index of the same name that referenced another feature class.

Although I am testing on a file geodatabase my target geodatabase for this workflow is in SQL Server and accessed via ArcSDE.


texts: Looking for real-world DJI Phantom 3 Professional sample geo-tagged images for testing

texts: I am looking for sample DJI Phantom 3 Professional imagery for testing purposes with the Drone2Map for ArcGIS Beta software. So far I am finding the application very polished and easy to use but I am running it on the sample datasets provided by ESRI which maybe be very clean and "ideally" suited for this software.  I would like to purchase the above-mentioned UAV do do 3D models and simple volumetric calculations for small sites. I would like to make sure that the software (which will no doubt be a very expensive yearly subscription) does in fact work as well with real world data. I would also like to run some test on some open source software packages some of which are mentioned here.


texts: Which projection should be chosen to model transmission lines across the USA?

texts: I am looking at different options for high voltage power lines in the US. These power lines cross several states. Since in the end the distance is of importance, is the USA Contiguous Equidistant Conic Projection the right projection? If not, which one would you suggest?


texts: basic problem using division in field calculator

texts:
This is a very basic question, but I can't figure out what's going on. I've got election results I'd like to simply display as percentages rather than vote totals, which should be as simple as dividing the column "ElectionYes" by "ElectionTotal". When I write that expression ("ElectionYes" / "ElectionTotal") into the field calculator or layer display box, though, it spits out "0" for all the rows! What gives? I should be seeing returns ranging from .402 to .670.

I've made sure that all the columns are integers, all the data is saved, and there are no null values in either column. All the other mathematical functions work (as in, I can create a column of "ElectionYes" * "ElectionTotal" with no errors), but I can't figure out how to do it with division. Any ideas?


texts: K-shortest path between two points in QGIS

texts:
I am new to QGIS. I have a road network layer in QGIS. I can compute the shortest path between the two points over the road network layer via the Network Analysis toolbox in QGIS.
I am interested in generating k-shortest paths between two points. I would like to ask if QGIS offers this functionality or is it possible to use additional Plugins/Python scripts to do the job.


texts: QGIS 2.12 Lines to Polygons not work properly

texts:
In QGIS 2.12.3 I have a small shapefile with lines (as often occurs in parcel boundary datasets) of 2 vertices each, snapped to each other (no topology errors, no dangles, no intersections). But the "Lines to Polygons" (screenshot) tool doesn't produce any polygons (in contrast to the two expected polygons).

What am I doing wrong or what am i having to do to make this tool work properly? (I'm checking out 2.14. tomorrow...) I know that in PostGIS the ST_Polygonize works pretty well, but importing my data to PostGIS first or any other workaround seems not a staightforward solution to me.


texts: How many polygons does a single polygon intersect with?

texts:
I currently possess a shapefile that consists of >5000 polygons. I would like to know, for each polygon, how many other polygons it intersects with and the names of each polygon it intersects with?

My output should be something like this:

Based on this previous question and @JeffreyEvans suggestion on using gIntersects, I have tried the following below. But, the process has been running continuously without stopping for the last 5 hours. Any suggestions?

Here's the code I ran in R:

require(rgdal)
library(maptools)
library(rgeos)
library(GISTools)
Poly <- readOGR("C:\\Users\\rameshv\\Terrestrial_dissolved_multipart.shp")
Poly_join <- gIntersects(Poly)

EDIT

I was able to solve part of my own question. However, I am still unclear as to write a list of comma separated names on each row. As of now, I have a loop that read in every shapefile, checks for an Intersection with a large Spatial Polygons Data Frame and prints out a dataframe of  and number of polygons it intersects with.

library(rgdal)
library(GISTools)
library(sp)
library(raster)

##List all the shapefiles
shapes <- list.files(path =    "C:\\Users\\rameshv\\Desktop\\Output_All\\", pattern=".shp$", full.names=TRUE)
head(shapes)

##Read in the big spatial polygons dataframe - this contains ALL MAMMALS
Big <-readOGR("C:\\Users\\rameshv\\Desktop\\Ranges_Overlap\\multipart.shp")

for (i in 1:length(shapes)){
    a <- readOGR(shapes[i])
    o <- over(a, Big,returnList = T)
    s <- o[[1]]
    count <- length(s$binomial)-1
    a$binomial <- as.character(a$binomial)
    results <- data.frame(a$binomial,count)
 write.table(results,file="C:\\Users\\rameshv\\Desktop\\overlapfile_Final.csv", row.names=F,append = TRUE, col.names = F, sep = ",")
}

Any thoughts?


texts: Determining which US zipcodes map to more than one state or more than one city?

texts:
I'm using a zip code listing, and I am curious to know how many (or which) zip codes map to more than one US state or US city?

For instance I know that Zip Code 42223 resolves to US Army, Fort Campbell which straddles the KY-TN state line. Oddly, the google API only returns TN for that state corresponding to that zip.


texts: What does the 'Field' selection do in QGIS report composer?

texts: In QGIS 3.4 report composer, when making a field group section, there is a selection under 'layer' labelled 'field'. Here, one seems to be able to select from all the fields in the layer. But as far as I can tell, this selection doesn't actually do anything. The report produces one map zoomed to each feature in the layer regardless of what field is selected here. Tables filtered by atlas id still produce one result per feature and there's no @layout_id variable. I can't fathom what the 'field' option is doing here.


texts: My organization uses several address locators stored in SDE.  The server license was recently updated, and shortly after, some of these locators simply disappeared -- they are no longer visible or accessible to ArcCatalog, Python, etc., and the underlying SQL Server shows no corresponding tables.  However, trying to create a new locator with the same name fails with the message that it "already exists."  Has anyone experienced a similar issue or have some other suggestions to either delete or restore these missing locators?


texts:
I'm trying to import a DWG file into QGIS. I tried with both 3.16 and 3.4 versions. Not able to open in either version. Attached is the screenshot of the error.

Downloaded a sample DWG from internet. Even the sample file also throwing the same error.


texts:
I often use closed lines to cut polygons on another layer.  When a line is not closed the cutting operation fails.  Is there an out-of-the-box tool that can easily tell me that a line is not closed and identify the vertices that need closing?

EDIT: Forgot to mention I have a standard license so the solution should work with a standard license.


texts:
I'm currently implementing 3DTiles into my java application. I know that I will only work with the region BoundingVolume. The BoundingVolume has 6 values: west, south, east, north, minHeight and maxHeight. Min and maxHeight are in meters above the WGS 84 ellipsoid and the other values are given in WGS 84 datum as defined in EPSG 4979 and are in radians. So far so good.
Now my problem is the next step. I am personally working with Point Clouds, so I read in a .pnts file and obtain the position data for every point that is in the Point Cloud. I get x,y and z values. I think the z values are in meters again but I don't know how to use the x and y values. If I try to convert them from radians to degree the values are not correct. I also get the values of an RTC_Center which has 3 values. I guess they are x,y and z but I also don't know how to work with these. I want to display them on a globe and need them as Latitude and Longitude(atleast I think so).
So my question is: how do I work with the obtained position values if I need them as Latitude and Longitude. I have the feeling the RTC_CENTER is important but I don't know for what.
I will leave some links to the 3DTiles specifications that I think will help.
Coordinate reference system
PointCloud point positions


texts: ArcObjects - Listen to Editor events in Extension Add-In using Java

texts:
I'm struggling at finding a way to listen to Editor events in my extension add-in. What am I trying to do is: after user creates a new point feature, extension should populate an attribute of it with its Z value. So I guess my extension needs to be aware of IEditEventsOnCreateFeatureEvent dispatch, but I don't know how to handle it using ArcObjects for Java. Documentation and forums are full of suggestions using VB and C#, and for ArcGIS Engine, but no Java Extension Add-In AFAIK.

Everything I've been able to do until now is getting an IEditor reference (which is not much, I know).


texts:
My application is currently displaying some search results in a dojo grid inside a dialog. I've been asked to embed a map inside the dialog but am struggling to do so. Has anyone done this?

I've got a simple map inside a jsp and have tried

dialog.set("href","url_for_my.jsp")

Also tried putting a button on the page to with map.render("mapdiv") after the dialog is displayed. Dialog.refresh() has no effect. There are no errors and no requests made for map tiles.

Can anyone suggest a better way of displaying a map in a popup?
Thanks!


texts: Sentinel 2 imagery does not align with Landsat imagery regardless of coordinate system

texts: I have downloaded Sentinel 2 imagery of Antarctica and imported it into a QGIS project that also contains older Landsat imagery, line vectors, and a basic basemap. The Sentinel 2 imagery was automatically downloaded in EPSG:3857 but the other layers are all in EPSG:3031. Regardless, the Sentinel 2 imagery appears in the correct location on the basemap. But in comparison to the Landsat imagery, the Sentinel imagery appears misaligned by roughly 120m. When the projection of the Sentinel data is changed to EPSG:3031 the imagery moves north to roughly the middle of the Pacific Ocean instead of becoming perfectly aligned with the Landsat imagery. Does anyone know why these two different sets of imagery in different coordinate systems are very closely, but not perfectly, aligned, and when put in the same coordinate system become very misaligned? If so, how do I fix this so both sets of data become aligned?


texts: Comparing data from different sources using QGIS?

texts:
I have two files in QGIS about restaurants in the same area.

How can I get elements that are only in the first source and not in the second and viceversa, connected to the same activity?

I also need to get the elements in common, associated to the same restaurant. The problem is that the attribute tables of the two sources are organised in different ways, there's not an univocal ID to join the sources, also the attribute "names" are slightly different between them and same cells are empty in the second source.


texts: Simple Spatial Join of Point Feature to Polygon Target Feature Not So Simple (Source Files Included)

texts:
I have what should be a simple task to perform: join ~47K census blocks from the PL redistricting files with pop and racial data joined to Chicago Public Schools High School boundaries to get an aggregate total of  >18 population.

I converted the blocks from polygon to point using centroids to avoid blocks from intersecting two boundary attributes in the target feature and thus being counted twice.  However, when I run the spatial join tool with Completely Contains, One_to_One aggregate, I still get an overcount of approximately 10k blocks and 700K people. I know there is a simple problem underlying this, so I'm hoping someone can spot the obvious and help me out here.

I've attached a dropbox link here:Direct Dropbox URL of Files to the two files that need to be joined for further clarification.

All Districts Merge refers to the target for the join, ...Centroids, the blocks to be joined to the Districts.

As an added bonus, if anyone knows of any FOSS redistricting tools--my end goal is to create 12-13 amalgamated districts composed of these boundaries.

Edit: I'm looking to summarize the pop total. That attribute is joined to the point file. It's a tabulation of pop over 18.

Before spatial join it amounts to approx 2,100,000, which is correct for Chicago city. There should be~47 k blocks too. When joined to the districts file however I get a sum of all districts of 2,800,000 and a join count for all districts of 57k.

This leads me to think that when joining blocks to districts I'm getting a double count of blocks that might slightly overlap two districts.

Even after setting join to completely contains parameters I had the same overcount.

Can you do a join where each feature is counted once and assigned to only one district?


texts: How to convert lat long coordinates into categorical value in R?

texts:
I have latitude, longitude coordinates as follows

latitude      longitude
41.775974     -71.329887
42.184913     -71.9179
41.682957373  -71.56037000

I want to use it as one of the features in a RandomForest classification. How to use it as a meaningful format?

Why i need this to be a standard format is Random forest is a classification algorithm if we use lat, long coordinates as a feature for classification means, the number or somewhat similar the difference is .6 or .5. so we cannot go and use lat long coordinates directly as a feature


texts: Trouble using saveACopy with group layers

texts:
This code works fine on layers in the TOC that are not group layers, but errrors out if there are any group layers in the TOC.

import arcpy, re, os
arcpy.env.overwriteOutput = True

#loop through layers and save a copy

def layer_save(mxd_path,path):

    mxd = arcpy.mapping.MapDocument(mxd_path)
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    fcs = arcpy.mapping.ListLayers(df)
    for item in fcs:
        formatted = str(item).replace("/","-")
        item.saveACopy(str(path) + "\\" + formatted + ".lyr")

if __name__ == '__main__':

    #variables
    mxd_p = arcpy.GetParameterAsText(0)
    path_lyrs = arcpy.GetParameterAsText(1)
    #kickoff
    layer_save(mxd_p,path_lyrs)


texts: Calculate the mean for all the attributes of a zipcode?

texts: I have a spreadsheet with two columns: Zipcode and Housing Price. Is it possible to calculate the average Housing Price per zipcode without manually using Excel commands per zip code?


texts: How can I export a shapefile from PostGIS with a GeoJSON query?

texts:
I tried to export a shapefile from my PostGIS installation using pgsql2shp as follow

pgsql2shp -f tiles.shp -h <host> -u <user> -P <pass> <database> "SELECT id, the_geom FROM <table> WHERE ST_GeomFromGeoJSON('{ "type": "Polygon", "coordinates": [ [ [ -92.472398018272358, 18.086381878379395 ], [ -92.472398018272358, 21.705890714546868 ], [ -86.679183238981381, 21.705890714546868 ], [ -86.679183238981381, 18.086381878379395 ], [ -92.472398018272358, 18.086381878379395 ] ] ] }') && the_geom;"

and the command gives the error

Error executing user query: ERROR:  quoted object property name expected (at offset 2)

I suppose is a problem related with the GeoJSON query because when I run

pgsql2shp -f tiles.shp -h <host> -u <user> -P <pass> <database> "SELECT id, the_geom FROM <table> WHERE id = 1;"

the shapefile is generated without error.


texts: Accounting for Spatial Autocorrelation in ZINB Model?

texts:
I'm doing an analysis on disease in Chile. I'm using the rate as the Dependent Variable and both remote sensing and housing-type data as the Independent Variables. I have several questions regarding my methods:

The data appears to follow the zero-inflated negative binomial distribution. I checked with the countfit tool in Stata, and the data best fits this distribution. Is it ok to fit a ZINB model to continuous data?
I need to account for spatial autocorrelation in both the dependent, and the remote sensing data. I'm not sure exactly how to do this. After some reading, I thought that I needed to generate a spatial weights matrix (using ArcMap). If this is correct, how do I implement it? I generated one for the disease rate, and used it in Stata: [pweight = weight], which greatly improved the model fit. However, I'm not sure if this is the correct usage.

Can anyone give me some pointers, or readings, or anything? I really want to complete this project!


texts: How do I group multiple polygons into the same record?

texts:
I am trying to figure out how to group multiple polygons into the same record (using something like a spatial join tool) using CartoDB.

I'd prefer to keep the polygons as separate rather than joining them to create one larger polygon, but I need them to be part of the same record


texts: Create random points, only where it overlaps specific values in a raster

texts:
I need to create a specific number of random points THAT FALL ON A RASTER with values between 0 and 1. The trick is the raster has a lot of holes - it is a grid of linear features where there are height values on the linear features, but everything in between is a NA. I have selected all my used locations that fall on the linear features in the raster, and intersected them to get the height value. Now I need to create the same number of random points, restricted to being where there are raster values (i.e no points are generated in the holes, where there are NA values).

Any advice? I've considered converting my raster to a polygon, however the raster is over a large extent (60 x 90 km) and is very high resolution (1m pixels), so this would probably take a long time.


texts: Installing QGIS 2.14.9-Essen with Qt 5.6.1 on RHEL7 (64-bit)


texts: No fetch Python plugins in qgis 1.8 + ubuntu 13.04

texts: I tried to install qgis 1.8 to ubuntu 13.04 but when installation finished I couldn't find "Fetch python plugin" button in plugin menu. I didn't find any soloution for this problem. Any idea what's the problem?


texts: Understanding Geodatabase Archiving and H tables in ArcSDE?

texts:
We have ArcSDE 10.0 on MSSQL 2008 with Geodatabase archiving enabled.  In the process of updating a topology I deleted the topology, disabled archiving, unregistered as versioned and then recreated the topology.

It seems that when I did this a copy of the feature classes in the feature dataset were created at the root of the geodatabase suffixed with an "_H".  When I deleted and created the topology again the features classes were copied to the root of the gdb again with an "_H1" suffixed.

Does anyone know what's going on?  I was told these were the Archiving tables.  Please provide links if available.


texts: Is there a Post-classification Majority Analysis feature on Google Earth Engine (Similar to ENVI)?

texts:
In ENVI there's a Majority/Minority Analysis function that allows you to to change spurious/"wrongly classified" pixels within a large single class to that class (it also allows you to enter kernel size and the weight of the center pixel in the kernel).

I have a classified image but I want to fix the "salt-and-pepper" effect of pixels the output image has.

Is there a Majority analysis function post-classification in Google Earth Engine?

What are sample scripts, if ever? And/or where would I start reading up on it? I can't seem to figure out where it is in the Earth Engine tutorial.

If there isn't are there any work-arounds to it? I'm trying to only use freely available software. Is there a function in QGIS perhaps?


texts: Display WMS Layer With Geoserver and Leaflet

texts:
I am trying to display a WMS layer from geoserver to my Vue.Js application. I am using Vuelayers (Leaflet) for the map rendering.

This is my code:

   <template>

WMS Layers with Layer selector

        NOTE: Use layer symbol on the top right to select which layer you
        want to display. More information about WMS (Web Map Service) can be
        found on the
        leaflet.js WMS example page.

    <l-map
      :zoom="zoom"
      :center="center"
      style="height: 100%"
    >
      <l-control-layers />
      <l-wms-tile-layer
        v-for="layer in layers"
        :key="layer.name"
        :base-url="baseUrl"
        :layers="layer.layers"
        :visible="layer.visible"
        :name="layer.name"
        layer-type="base"
      />
    </l-map>

</template>

I don't have any error codes and nothing it's displaying in my browser. The integration is bit different in vue.js and I'm a bit confused. Any clues where might be my mistake?

My link layer: http://localhost:8081/geoserver/cite/wms?service=WMS&version=1.1.0&request=GetMap&layers=cite%3Avnm_polbnda_adm3_2014_pdc&bbox=102.144027709961%2C8.4089241027832%2C109.471290588379%2C23.3922080993652&width=375&height=768&srs=EPSG%3A4326&format=application/openlayers

This is the my developer screenshot :


texts: Access a topology sample code runtime error when define ITopologyExtension

texts:
I encountered a runtime error when I used the sample code from ESRI on Accessing a topology. The sample code is as follows:

Dim topoUiD As UID = New UIDClass
topoUiD.Value = "esriEditorExt.TopologyExtension"
Dim m_application As IApplication   ' this line is added by me
Dim topologyExt As ITopologyExtension = CType(m_application.FindExtensionByCLSID(topoUiD),  ITopologyExtension)
Dim topology As ITopology
topology = CType(topologyExt.CurrentTopology, ITopology)

A runtime error occurs when execute:  Dim topologyExt As ITopologyExtension = CType(m_application.FindExtensionByCLSID(topoUiD),  ITopologyExtension

I guess it is caused by the definition of m_application, anybody can tell me how to define it? thanks.


texts: Selecting and updating from multiple tables using arcpy.UpdateCursor()?

texts:
I have two feature classes which have common Valve ID #s.  One table has multiple service connections with a unique number.  I need to select all services linked to each valve and collect the individual service # then update the service fields (1-4) in the Valve class with each connection #.

I believe I can do it with an updateCursor but I'm not a strong coder.


texts: How are GeoTIFF pyramids/overviews standardised?

texts:
While pyramids/overviews are not part of the GeoTIFF standard, many tools support creating them. For example vips/nip2, Orfeo Toolbox (otb) and ossim which all promise some support for creating those. But I have no idea if they would all create files in a format that the others support. Looking through the documentation of the tools does not really mention anything about that.

vips is not a geospatial tool and I could not find any enduser-friendly documentation but the IIPImage docs mention its support for "Tiled Pyramidal TIFF" : http://iipimage.sourceforge.net/documentation/images/

otb does not mention anything about a format or spec for its "Multi Resolution Pyramids": https://www.orfeo-toolbox.org/CookBook/CookBooksu65.html

ossim says it supports various output formats for its "reduced resolution data sets" but I don't know what they actually mean: https://trac.osgeo.org/ossim/wiki/img2rr

GDAL also does not really specify things about its "overview images": http://www.gdal.org/gdaladdo.html

So they all have pyramids/overviews but it is not clear if they are cross-compatible.

On more general pages I found the following quotes:

http://iipimage.sourceforge.net/documentation/images/ says

  Tiled Multi-Resolution (or Tiled Pyramidal) TIFF is simply a tiled
  multi-page TIFF image, with each resolution stored as a separate layer
  within the TIFF. This is a standard TIFF extension and is supported by
  most image processing applications including Photoshop, GIMP, VIPS and
  ImageMagick. The libtiff codec library is also perfectly capable of
  reading and writing such images.

Is that what everyone uses?

The Library of Congress also has some info: http://www.digitalpreservation.gov/formats/fdd/fdd000237.shtml
They note:

  Pyramid TIFF files created by different applications are not
  necessarily the same in structure. In particular, judging from
  analysis with JHOVE and the identify command in ImageMagick, Adobe's
  Photoshop and Image Magick generate files with different internal TIFF
  structures; in both cases, most software that can handle TIFFs appears
  to recognize the primary (full-size) TIFF without problem.

So, are those formats standardised and specified and documented somewhere? How can I find out which tool can produce them in compatible ways to others? Are there any geospatial tags in overviews/pyramids or could I use any software to create them on the image data?


texts: ogr2ogr - Looping through files in a folder on a Mac

texts:
I am attempting to use OGR2ORG on Mac OS X to loop through 30 shapefiles, change their projection and save the files into a new directory. I am new to using Macs and the OGR2OGR command line environment.

The original files are in /data/ChicagoBaseLayers/clipped.  I would like the reprojected layers to go to /data/ChicagoBaseLayers/Projected.

After navigating to the 'clipped' directory on the command line, I tried

for %f in (*.shp) do ogr2ogr -t_srs EPSG:3857 -s_srs EPSG:4326 /Users/stevenchong/data/ChicagoBaseLayers/Projected/%f Users/stevenchong/data/ChicagoBaseLayers/clipped/%f

but I get the error message: syntax error near unexpected token `('

What am I doing wrong?


texts: How to select overlapping/non-overlapping polygons from a single shapefile in QGIS?

texts:
I have a shapefile with a series of polygons, some of which intersect each other and some that don't. I want to select only those that do not overlap with any other polygon.

I've tried the Select by Location tool with Intercept but this selects all polygons in the shapefile.

As an example, I'd like to select the two polygons to the right, since they intersect, but not the one to the left. Of course, I can't do this by hand since there are too many polygons.

I'm working on QGIS 2.10.


texts: ArcGIS 10.1 arcpy Update Cursor is returning a list instead of a row

texts:
I keep getting the following error:

  Traceback (most recent call last):
  File "R:\GIS_ScriptsandTools\PullHotParametersLabel.py", line 16, in

      TitleLine = ""+row.getValue(IDfield)+" ("+ row.getValue(Datefield)+")
  "

  AttributeError: 'list' object has no attribute 'getValue'

When I am trying to run a typical for row in cursor loop when working with an UpdateCursor.
For the life of me, I can't figure out why it thinks the row is a list and not a row object.
Here is my script:

import arcpy
SourceTable = arcpy.GetParameterAsText(0)
IDfield = arcpy.GetParameterAsText(1)
Datefield =arcpy.GetParameterAsText(2)
SelectedParams = arcpy.GetParameterAsText(3)
OutputField = arcpy.GetParameterAsText(4)
SelectedFields = SelectedParams.split(';')
AllFields = arcpy.ListFields(SourceTable)

FieldDict = {}
for field in AllFields:
    FieldDict[field.name]=field.aliasName

with arcpy.da.UpdateCursor(SourceTable,"*") as cursor:
    for row in cursor:
        TitleLine = "<und><bol>"+row.getValue(IDfield)+" </bol>("+ row.getValue(Datefield)+")</und>"
        LabelList = [TitleLine]
        for Field in SelectedFields:
            if row.getvalue(Field) <>0:
                Alias = FieldDict.get(Field)
                ParamValue = row.getValue(Field)
                ParamLine = '{0} = {1:,}'.format(Alias,ParamValue)
                LabelList.append(ParamLine)
        if len(LabelList)<2:
            LabelList.append("All &lt; 2L")
        Label = '\n'.join(LabelList)
        row[OutputField] = Label
        cursor.updateRow(row)
del row, cursor

Does anybody have any ideas where I'm screwing up on this?


texts: Understanding if-else statements in ArcGIS field calculator when coded value domain used?

texts:
I am attempting to set a newly created field based on the strings in another field 'TNMFRC'.  I am pretty sure that I have the code write, and yet everything is set to a value of '6'.  Here's my code:

 def myCalc(road):
  if road == "Local Road":
    return 1
  elif road == "Controlled-access Highway":
    return 2
  elif road == "Secondary Highway or Major Connecting Road":
    return 3
  elif road == "Ramp":
    return 4
  elif road == "4WD":
    return 5
  else:
    return 6

and in the other box:

myCalc(!TNMFRC!)

Here are screenshots of the fields and the field calculator:

Can anyone see what I might be doing wrong?  I can't figure it out.

EDIT:  So even though the field appears to be as string based field, it seems in the properties to be listed as a long integer numeric field?


texts: Download data from Geomet using WCS request

texts:
I'm trying to download some temporal weather raster data from a government supported server called GeoMet. I'm able to download some data using the request below:

http://geo.weather.gc.ca/geomet-beta?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&COVERAGE=GDPS.ETA_NT&BBOX=(-160,40,-30,80)&CRS=EPSG:4326&WIDTH=560&HEIGHT=900&FORMAT=GeoTIFF_16

But I still a bit confused about the request structure, and how to acquire them? ( for example, how to find out which version to use, how to find out the coverage name etc)．


texts:
I wonder if it is possible to check the geometry of multiple polygons using Fiona or shapely library In Python. Whichever polygon does have proper geometry will be kept and the rest will be deleted and a list of polygon names with whether they have correct geometry and not correct geometry be printed.
Problem Statement:
I have more than 200 polygons, out of these polygons, some polygons don't have the correct geometry, i.e. in the below-attached image ID, 19200102 polygon doesn't have geometry information.  I want to filter out, want to keep only which polygons have the correct geometry. Print the name of polygons correct and not correct geometry.


texts: arcpy Union_analysis general function failure problem

texts:
As part of a longer script, I use this code:

import arcpy, os

Polygon = r"pathtoashapefile.shp"
workspace = r"path_to_a_folder"
output = os.path.join(workspace, "test.shp")
arcpy.Union_analysis(Polygon, output, "ALL", "#", "GAPS")

and I get the error:

Runtime error
Traceback (most recent call last):
  File "test.py", line 6, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\analysis.py", line 563, in Union
    raise e
RuntimeError: Object: Error in executing tool

And this is a super frustrating error because it gives me no additional information.  I am only unioning with one polygon shapefile because I want to see overlaps between polygons within the shapefile.  I have tried all of the following:

putting Polygon in a list like: [Polygon]

Putting Polygon in a list like: [Polygon, "#"]

putting Polygon in a string formatted: string = '"%s" #' % Polygon

I have checked that Polygon has data in it and is not null.  I have checked that the path to Polygon and to workspace are valid.  I have not gotten any error other than the 'line 563' error in all of my attempts.


texts: QGIS: tiles from WMS server rendered incorrectly (maybe incorrectly made XML config?)

texts:
I want to load some tiled map in QGIS. The tiles are okay (I've checked them with leaflet map viewer), but when I load them into QGIS (adding WMS config XML file as raster layer), the tiles are rendered weirdly.

Here's my XML config file (chelyabinsk.xml):

<GDAL_WMS>
 <Service name="TMS">
  <ServerUrl>file://path_to_tile_folder_chelyabinsk2_/${z}/${x}/${y}.png</ServerUrl>
  <SRS>3857</SRS>
 </Service>
 <DataWindow> <!--
   <UpperLeftX>-20037508.34</UpperLeftX>
   <UpperLeftY>20037508.34</UpperLeftY>
   <LowerRightX>20037508.34</LowerRightX>
   <LowerRightY>-20037508.34</LowerRightY> -->
   <UpperLeftX>6621281.138</UpperLeftX>
   <UpperLeftY>7589891.161</UpperLeftY>
   <LowerRightX>6816959.931</LowerRightX>
   <LowerRightY>7301264.942</LowerRightY>
  <TileLevel>14</TileLevel>
  <TileCountX>1</TileCountX>
  <TileCountY>1</TileCountY>
  <YOrigin>bottom</YOrigin>
 </DataWindow>
 <Projection>EPSG:3857</Projection>
 <BlockSizeX>256</BlockSizeX>
 <BlockSizeY>256</BlockSizeY>
 <BandsCount>4</BandsCount>
 <Cache/>
</GDAL_WMS>

And here are the tiles: https://www.sendspace.com/file/vvvytc

Can you tell me, how to fix my XML?


texts: Python and Google Maps API place search: Read csv read row by row and loop fuction

texts:
I want to automate a task in Python with the Google API Search Place. Right now, my code with Google Maps API is returning places 'hospital' for one set of geographic coordinates.

But, I have a CSV file which contains many geographic coordinates.

1: Hanoi 10.762622, 106.660172
2: Ho Chi Minh 12.762622, 108.660175
3: Ho Chi Minh 11.8542, 108.660175
4: ...
5: ...

As you can see in my code, this is not efficient because I need to change the geographic coordinates manually.

For each row (geographic coordinates) in my CSV file, I want my code to read the geographic coordinates 1: Hanoi and give me results 'hospitals' for this geographic coordinates. After read row 2: Ho Chi Minh and same for all the other rows.

How I can achieve that or any good examples of making tasks like this more efficient?

import urllib
import urllib.request
import json

googleGeocodeUrl = 'https://maps.googleapis.com/maps/api/place/textsearch/json?query='
keyword = "hospital"
geolocation = "&location=-12.135,-77.023&radius=5000"
APIKEY = '&key='+'apikey'

url = googleGeocodeUrl + keyword + geolocation + APIKEY
print(url)

url = googleGeocodeUrl + keyword + geolocation + APIKEY
json_response = urllib.request.urlopen(url)
search = json_response.read().decode('utf-8')
searchjson = json.loads(search)

export = open('hopital.csv','w')
for place in searchjson['results']:
    print(place['name'])
    print(place['geometry']['location'])
export.write(place['name']+','+str(place['geometry']['location']['lng'])\
 +','+str(place['geometry']['location']['lat'])+'\n')
export.close()


texts:
When i try to do mutliprocessing in python Qgis it doesn't work

from multiprocessing import Process,Lock
import multiprocessing

def name1(l,i):
    print i

if __name__=='__main__':
    lock=Lock()
    for i in range(5):
        print (i,'process start')
        process=multiprocessing.Process(target=name1,args=(lock,i,))
        process.start()

The above code it's working perfectly in eclipse and other IDE's but it is not working with qgis python script

output from eclipse:
(0, 'process start')
(1, 'process start')
(2, 'process start')
(3, 'process start')
(4, 'process start')
1
2
3
4
0

But in pyqgis there's no error and no output


texts: How to calculate distance to a specific slope value (in degrees) from each data point on QGIS?

texts:
I'm using QGIS. Before this, I managed to get an attribute table showing depth/slope of each data point. Now, I’d like to add an additional column that shows the distance of each point to a slope value. So say for example I choose to specify a slope value of 1 degree, what would the distance be from point A/B/C etc to the nearest slope that is 1 degree?

But first I need to calculate the distance to slope from each data point which I can't seem to figure out. My slope values were derived from terrain analysis.


texts:
OSM Copyright and License dictates that the users will give OSM credit.

Is there a way to add the credit into QGIS maps automatically?


texts: QGIS Tool Scripting with "Select by attribute"?

texts:
Why am I getting the error message

  The following layers were not correctly generated:  Output

when writing a script that uses the QGIS geoalgorithm "Select by attribute"?

Here's my script (My "Input" and "Output" parameters are set in the GUI when the script is run.):

processing.runalg('qgis:selectbyattribute', Input, "STRNAME", 0, 'Robie', Output)

I used this format after looking at the HELP tab of the "Select by attribute" tool and seeing the following instruction:

processing.runalg('qgis:selectbyattribute', input, field, operator, value, output)

Confusing then, to find the answer (helped by @StephenKay) is to remove the final parameter ('output').  Would have been less misleading if the tool's own HELP tab did not show it as the 5th and final parameter, or at the very least indicated that it is an optional parameter (as the tool clearly works without it).


texts: What is the max amount of connections a GeoPackage can handle in GeoServer?

texts:
I'm currently working on a GeoServer deployment and for simplicity, I want to to transition to using GeoPackages rather than Shapefiles (longer column names, less work integrating the data, etc.).

As GeoPackages are based on the SQLite3 standard there's a possibility that, with too many connections, users won't be able to read the polygons saved within a GeoPackage. (Here's a similar discussion on StackOverflow).

I'm not sure whether GeoServer loads shapefiles and other data into a database that can accommodate 100s of users, or if it simply reads shapefiles/other vector data directly when a user demands access to them. While other SE's may know more about DBMSes, I felt this question was appropriate for the GIS SE users who have more likely used GeoServers.

My question is what are the connection limitations for GeoPackages served through GeoServer?


texts: Loading WMS layer from USGS WMS server using Leaflet?

texts:
I am trying to load WMS layer using leaflet from USGS WMS server but I am getting an empty image. I have gone through Loading WMS layer from external server in Leaflet? link but it doesn't solve my issue.

My code:

var map = L.map('map').setView([47.608013,-122.335167], 5);
    L.tileLayer.wms(USGS server link,
 {
  layers: '0',
  format: 'image/png',
  version: '1.3.0',
  transparent: true,
  continuousWorld: true,
  attribution: "",
  tiled:true
 }).addTo(map);

I even tried changing version to 1.1.1, removed continousWorld: true etc.(tried almost everything).  latitude and longitude that I used is Seattle's. Other examples of WMS Leaflet are working for me. The URL it is using to get the image is:

<img src="<wms server>?request=GetCapabilities&amp;service=WMS&amp;SERVICE=WMS&amp;REQUEST=GetMap&amp;VERSION=1.3.0&amp;LAYERS=0&amp;STYLES=&amp;FORMAT=image%2Fpng&amp;TRANSPARENT=true&amp;HEIGHT=256&amp;WIDTH=256&amp;CONTINUOUSWORLD=true&amp;TILED=true&amp;CRS=EPSG%3A3857&amp;BBOX=-13775786.985667603,6261721.35712164,-12523442.714243276,7514065.628545966" style="height: 256px; width: 256px; left: 270px; top: 17px;">

The USGS WMS server is working fine in QGIS tool.


texts: Reprojection of Landsat 8 using the existing shapefile in ArcGis

texts:
I have downloaded a landsat8 and when I overlay it with a shapefile of the location where I want to extract by mask , the combination of landsat 8 and the shapefile are not matching. The shapefile is outside of the landsat image. I reprojected the landsat8 using the same projection system but the problem is still there.

Is there another way I can fix that problem?


texts: Automatic Land Cover Classification Method in eCognition?

texts:
I'm exploring several methods of classifying multi spectral data in eCognition.

The classification will be fairly high level; I'm not trying to find individual species but want to identify areas of forestry, urban areas and water bodies.

Most methods require training areas, but in this instance, I need an automated method. eCognition appears to be the best solution as you can create rules based on the spectral characteristics of the land cover types but this is an expensive solution!

Has anyone used any alternative methods?


texts: Applying different labeling options at different scales to different points (that aren't easily separated from the other points)

texts:
I am creating maps (in ArcMap 10.1) for a web-based application.  There are more than 500 places to be labeled at some scale (the map scales range from 1:32 million to 1:1.5 million).  As you can imagine, some of the places are very near one another and some are far away from other places.

I would like those places that are farther away from other places (where there is no labeling conflict) to be labeled when the map is at a smaller scale than those for which zooming in to a larger scale is needed to avoid a cloud of overlapping names. (I don't mind them overlapping at the largest scales...)

I am using the Maplex labeling but don't see how to easily accomplish this. At some point they may still overlap, so I cannot use the "never overlap" option...

Thank you for any suggestions or guidance you can offer!.

Sandy


texts:
I was trying to convert the raster to vector/matrix using getValues and I am getting this error:
Error in (function (classes, fdef, mtable)  : unable to find an inherited method for function ‘getValues’ for signature ‘"character", "missing", "missing"

I have put the code below for reference.
#mosaic

Mosaic<- "Mosaic.tif"

#shapefile

shp <- "Tahsil Boundary2001.shp"

#Extacted Raster

gondia<- "extract.tif"

#convert raster to vector\matrix

nr <- getValues(gondia)


texts: Make a Declination Arrow for Data Driven Pages?

texts: I am making a map book over a large area and want to add a correct declination arrow for every page.  I have not found a way in ArcGIS to make a declination arrow (in 10.3).  I could make my own using two north arrows but I don't want to have to change it for everymap along the book. Is there a way to do this?


texts: Data for demonstrating world poverty on Mollweide projection?

texts:
I am overwhelmingly failing at finding data about global poverty that is any more specific than at the country level.

Because poverty has such a wide definition I am open to really anything that is an indicator, like high infant mortality or others.

Any suggestions to point me in the right direction for data more specific than the country level would be great?

Preferably it is not raster data but data I can manipulate on Arc 10.1.


texts: Sorting layers by name or type in QGIS Layers Panel

texts:
Oddly I'm not able to find an answer to how to simply sort layers by names or by type in the "Layers Panel". Can someone point me in the right direction?
QGIS version 3.14.15-Pi


texts: Making python script results appear in ModelBuilder?

texts:
I have a model (ModelBuilder) and within I added a python script to calculate a field. My problem now is, that I don't know how to make the output appear in my model as it is the same as the input (only with the field calculated).

I imagine there is a very simple way, but I am new to ModelBuilder and python and just can't think of it?

Alternatively, I could also try using the field calculator, but I would prefer it this way.

And here is my script:

import arcpy
from arcpy import env

env.workspace = arcpy.GetParametersAsText(0)  # the workspace
inLayer = arcpy.GetParameterAsText(1) # The shapefile
inField = arcpy.GetParameterAsText(2)  # The column in attr table with Area
outField = arcpy.GetParameterAsText(3) # The column in attr table which should be calculated

# Create a list from area field and find min & max value

listing = []
rows = arcpy.SearchCursor(inLayer)
for row in rows:
    listing.append(row.getValue(inField))
del rows

listing.sort()
min = listing[0]
max = listing[-1]

# Calculate the normalized value

rows2 = arcpy.UpdateCursor(inLayer)
for row2 in rows2:
    row2.setValue(outField, (max - row2.getValue(inField))/(max - min))
    rows2.updateRow(row2)
del rows2


texts: Sentinel-1A and Lansat-8 Resolution difference

texts: Sentinel-1A level-1 high resolution GRD data have 20X22m resolution (rgXaz)with 5X1 number of looks. While Landsat-8 data have 30m resolution. My question is how can I compare the resolution differences between them? In other words, which one of them has higher resolution? When I digitize a lake size from the two sensors, lake size  from  Landsat-8 is always greater than Sentine-1A. Note that its a same lake; images were taken on the same acquisition date or with not much of time difference. What could be the reasons behind for differences in the lake sizes? Is it because of resolution difference? If yes, How can I  interpret the results?


texts:
I am using Oracle to store ArcGIS data. I have created my table named PARCELS from ArcCatalog, and I inserted some data using ArcMap. So when I run the following SQL query: select * from MDSYS.ALL_SDO_GEOM_METADATA;

OWNER   TABLE_NAME   COLUMN_NAME   DIMINFO             SRID
MY_SCH  PARCELS      SHAPE         {{null, 12345.4,... NULL

In this stuation, I can not run transform queries. Why this issue is occurred? How can I fix this?


texts: Error while changing the symbol of a layer in QGIS in plugin

texts:
I am trying to add a feature then trying to change shape and size of feature like below:-

layer = QgsVectorLayer('Point?crs=EPSG:4326', 'point' , "memory")
pr = layer.dataProvider()
pr.addAttributes([QgsField("ID", QVariant.Int)])
pt = QgsFeature()
layer.startEditing()
point1 = QgsPoint(float(lon),float(lat))
pt.setGeometry(QgsGeometry.fromPoint(point1))
pt.setAttributes([inputid])
pr.addFeatures([pt])
layer.updateExtents()
layer.commitChanges()
QgsMapLayerRegistry.instance().addMapLayers([layer])
layer = qgis.utils.iface.activeLayer()
symbol_layer = QgsSimpleMarkerSymbolLayerV2()
symbol_layer.setSize(5.0)
symbol_layer.setColor(QColor(255,0,0,255))
layer.rendererV2().symbols()[0].changeSymbolLayer(0, symbol_layer)

but this give me error :-

layer.rendererV2().symbols()[0].changeSymbolLayer(0, symbol_layer)
TypeError: expected string or buffer


texts: What open source solution would work best for auto extracting features from satellite images?

texts:
Possible Duplicate:
  Feature Extraction - Land cover

I have been looking for an open source solution which would allow me to define a region (or collection of images) and extract features from such images.

I realize such capabilities are far from perfect however my hope is that there is some type of capability which exists by which I can manually define example features I want to extract and through some type of machine learning procedures they application can identify others which are similar based upon the color variations of the pixels.

My goal is to end up with a vector shapefile which the results of this extraction.

Any assistance or guided procedure would be very helpful here. I would really like to find some solution which works in conjunction with QGIS, GRASS or R.


texts: How to output 3d geometries via grass plugins v.drape, v.to.3d in qgis?

texts:
When I try to process either a 2d Linestring or a 3d MultilinestringZ vector through the v.drape or v.to.3d processing tools in qgis, I get an error stating that it is outputting a 2d geometry instead of 3d and should set a lco flag to force 3d geometry:
WARNING: Vector map <output660b8ef0560648d183efaf4701ffc6ea> is 3D. Use format specific layer creation options (parameter 'lco') to export <in 3D rather than 2D (default).

This checks out because the file (that is created successfully) reports Geometry: Line(Linestring) in the layer properties
So, what flag is needed to force a 3d geometry with v.drape? I looked at the v.drape docs but could not determine what the lco flag should be.


texts: Area Value not getting reflected in attribute Table using PyQGIS?

texts:
I have a WardBoundary.shp file in my project. I need to find the area of each polygon ( representing different wards).With the help of this comment on Adding Area column in Polygon layer using PyQGIS?, I executed the code. I didn't get any error but still my attribute table doesn't show area value for any ward. Only the column Area is added in the field list. Here is the link of my shape file

 #!/usr/bin/python3

  from qgis.core import QgsVectorLayer,QgsVectorDataProvider
  from qgis.PyQt.QtCore import QVariant

  # ****Choose**** where you want to calculate your area

  # from a layer stored somewhere in your computer
  layer = QgsVectorLayer('C:/Users/Prachi/Desktop/Qgis   sample data\Ward_Saphale_updated.shp', 'WardBoundary', 'ogr')

  # in the active layer in the TreeLayer (the undelying layer)
  #layer = iface.activeLayer()

  # Here we get the capabilities of your layer (Add attribute layer, edit feature ect ..
  caps = layer.dataProvider().capabilities()

 # We make a list of fields from their name
 fields_name = [f.name() for f in layer.fields()]

# We check if we can add an attribute to the layer.
if caps & QgsVectorDataProvider.AddAttributes:
# We check if the attribute field is not exist
if "Area" not in fields_name:
    # We add the field name Area and with the double type (it can be integer or text
    layer.dataProvider().addAttributes([QgsField("Area", QVariant.Double)])
    # We update layer's field otherwise we'll not have the field
    layer.updateFields()
    # Recreate the list field by the name to have index of the field
    fields_name = [f.name() for f in layer.fields()]
    # we get the index of the Area field
    fareaidx = fields_name.index('Area')
else:
    # We are here because there is a field name Area
    print("The Area field is already added")
    # Recreate the list field by the name to have index of the field
    fields_name = [f.name() for f in layer.fields()]
    # we get the index of the Area field
    fareaidx = fields_name.index('Area')

 # Here we check if we can change attribute of the layer
 if caps & QgsVectorDataProvider.ChangeAttributeValues:
  # we loop^on every feature
    for feature in layer.getFeatures():
    # For each feature :
    # We calculate the area and put the index of the field Area
    # We round the area value by 2 digit
    attrs = {fareaidx : round(feature.geometry().area(), 2)}
    # We change the the value of Area Field for this feature.
    layer.dataProvider().changeAttributeValues({feature.id() : attrs})


texts: Creating new GeoServer user through rest api?

texts:
I've tried posting this xml:

<user>
    <userName>Jim</userName>
    <password>password</password>
    <enabled>true</enabled>
</user>

to these addresses:

http://localhost:8080/geoserver/rest/usergroup/user
404 - "Not Found"
http://localhost:8080/geoserver/rest/usergroup/users/user
404 - "Not Found"
http://localhost:8080/geoserver/rest/security/usergroup/user
404 - "Not Found"
http://localhost:8080/geoserver/rest/security/usergroup/users/
500 - "Server Error"
http://localhost:8080/geoserver/rest/security/usergroup/users/user
404 - "Not Found"
http://localhost:8080/geoserver/rest/security/usergroup/user/Jim
500 - "Server Error"
http://localhost:8080/geoserver/rest/security/usergroup/users/Jim
404 - "Not Found"
http://localhost:8080/geoserver/rest/security/usergroup/users/user/Jim
404 - "Not Found"
http://localhost:8080/geoserver/rest/security/usergroup/service/MyGroup/user
404 - "Not Found"
http://localhost:8080/geoserver/rest/security/usergroup/service/MyGroup/users
500 - "Server Error"
http://localhost:8080/geoserver/rest/security/usergroup/service/MyGroup/users/user
404 - "Not Found"
http://localhost:8080/geoserver/rest/security/usergroup/service/MyGroup/user/Jim
500 - "Server Error"
http://localhost:8080/geoserver/rest/security/usergroup/service/MyGroup/users/Jim
404 - "Not Found"
http://localhost:8080/geoserver/rest/security/usergroup/service/MyGroup/users/user/Jim
404 - "Not Found"

I've tried with and without the GeoFence plugin.

Is this actually possible?

I have succeeded in creating roles, usergroups, workspaces, coverageStores and raster layers from GeoTiff files, and dataStores and vector layers from shapefiles, but no progress on creating users.

I gather that some modules require you to change the admin password from the default, so I've done that too.  My admin account has ADMIN, GROUP_ADMIN, and ROLE_ADMIN roles, and I've even added a USER_ADMIN role in the case such a thing exists.

I'm running geoserver-2.12-SNAPSHOT from 18 September.


texts: Calculate total area of both adjacent and contiguous raster cells within a radial distance around a polygon's centroid

texts:
I have a raster containing grid-cell values of either 1 or no data, NA, representing used- or unused-habitat types, respectively.
Using R's raster package, I want to calculate the total area (in hectares) of a contiguous block of adjacent cells that fall within a defined (either square or circular) radial distance around a polygon's centroid.
My raster has a resolution of 30-m by 30-m (e.g., with a Lambert azimuthal equal-area projection), and each grid cell is approximately 90 square-meters in size. After counting the number of contiguous cells within a block, as defined by the radial distance around a centroid, I then multiply the total cell count by 90 (i.e. ignoring NA values) and later divide that product by 10,000 to convert from the total area in square-meters to the total area in hectares.
The problem that I am running into is that my usage of raster::ncell and raster::clump together returns the (grand total) number cells within a block, despite asking for no data, NA, values to be ignored. I also wanted the code to consider only the contiguous block of cells that is immediately adjacent to each polygon centroid within the defined radial distance. In other words, I did not want to quantify the total area based on other contiguous blocks separated from the immediately adjacent block by  the no data, NA, features.
Here's what I have tried so far:
cropped_locale <- crop(x=raster, y=extent(subset(centroids, Locale=="Name"))+5000) ##Cropping raster 5000 m around one centroid
total_area <- (ncell(clump(cropped_locale, directions=8, na.rm=TRUE))*90)/10000 ##Converting from square-meters to hectares

And here's an image of a subset of the raster data, including one polygon's centroid:

Perhaps my usage of both the ncell() and clump() functions is an incorrect approach to such calculations, and I am open to any suggestions.


texts: How to check ArcGIS REST Services Directory for date and time of publication?

texts: Is there a way to check quickly within the service description when (date, time) the service had been published?


texts:
Is Leaflet 100% free?

if you  go to doc example , you will see :

var mymap = L.map('mapid').setView([51.505, -0.09], 13);

L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token=pk.eyJ1IjoibWFwYm94IiwiYSI6ImNpejY4NXVycTA2emYycXBndHRqcmZ3N3gifQ.rJcFIG214AriISLbB6B5aw', {
    maxZoom: 18,
    attribution: 'Map data &copy; OpenStreetMap contributors, ' +
        'CC-BY-SA, ' +
        'Imagery © Mapbox',
    id: 'mapbox.streets'
}).addTo(mymap);

Which means that the tiles are supplied by mapbox.


texts: How to configure Orfeo toolbox in QGIS 2.18 on Windows 7 32-Bit?

texts:
I had installed QGIS 2.18.3 using OSGeo4W and then the OTB standalone binary.
In QGIS, I had manually entered the directory of the

OTB applications folder (C:\Program Files\QGIS 2.18\apps\OTB-5.10.0-win32\lib\otb\applications) and
Command line folder (C:\Program Files\QGIS 2.18\apps\OTB-5.10.0-win32\bin)

and I had checked the "Activate" button.

But nothing happens when I close the Processing window and there is no OTB menu appearing in the Processing Toolbox.

How can I resolve this issue?


texts: Creating fishnet with row Iterator and fishnet extent does not update to reflect each row feature

texts:
I'm a writing about a geoprocess that I can do on an individual basis but for some confounding reason am unable to do with the row iterator in Model Builder. Given that I want to repeat the same process 68 times, I hope that I can be spared the pain of doing this the long way (one-by-one). I have tried to be as complete as possible in describing the problem, including pictures and detailed commentary, hence the length of this post.

The task in question is exceedingly simple. Let's say there is one square. I want to create a fishnet bounded by the extent of the square. In the create fishnet tool, in the "template extent" drop down box, I select the file with the square. Each square happens to be 2200 meters by 2200 meters. So I create 22 rows and 22 columns of 100 by 100 meter cells. The entire square is covered by the fishnet.  So far so good. I can do this part.

Now let's say I have a shapefile that contains 68 rows. Each of these rows corresponds to one of 68 squares scattered across a metropolitan area (all the square are the same size). This is what I am talking about: one shapefile in the table of contents and 68 squares:

You will see that the attribute table is quite normal:

I want to do exactly what I have described above with one square, except now the file in question contains 68 rows instead of one and I want to use the row iterator function to execute the task on each row of the file.

I set up the process like so in Model Builder: the file with the 68 rows is the "input table" to the "iterate row selection" tool. If I understand correctly, each row is treated as a unique item by the darker green "I_FeatureEnvelope_68HHs." I subsequently want to treat each of these unique features as inputs to the "create fishnet" tool.

When I go to connect the "I_FeatureEnvelope_68HHs" to the "Create Fishnet" box I see a pop up box that says Template Extent. I click Template Extent, thinking that it means: "treat each feature represented by each of the 68 rows as unique extents to be used in creating 68 different fishnets." (sort of like if you were to open the attribute table, highlight a row and click 'zoom to selected' each row corresponds to an extent)

When I double click on the Create Fishnet box, I see this:

Template Extent has already been chosen for me. The selection appears to make sense. I set the cell width and height and numbers of rows and columns: 100m x 100m  and 22 x 22 as I stated before. I select polygon geometry type at the bottom.

I run this:

Let me tell you what I think I should see. To recap, there are 68 rows. Each of these corresponds to a square. I am creating 22 row x 22 column fishnets for each square which leads to a 484 cells fishnet per square. That means the output file should have 68 x 484, or 32,912 rows (or 484 cells per each square).

Let me tell you what actually happened. My process led to the creation of only one fishnet in the lower left hand corner:

Close up it looks alright, a 22 x 22 grid with 484 cells:

I have a faint idea of what is going on. Consider the file with the 68 squares. The one fishnet that was created corresponds to the lower left hand corner of the full extent of this file - so the Template Extent is not considering each of the 68 rows as unique items. It simply takes the full extent of the 68 row file and focuses on the lower left hand corner.

I have played around changing the settings a bit, but to document all that I'll end up writing a novel. To make a long story short, I have found it impossible to use the row iterate tool to create a fishnet that corresponds to the extent of each of the 68 features in the file in question.


texts: How to compress raster image for publishing?

texts:
I have a (8-bit) satellite image that is e 6 gb. Is there a way I can compress the image smaller size without losing the quality?

Format: .tif
Reason: After publish the image, I notice a higher usage or ram, so I assume it has something to do with the image.


texts:
I am looking to work with hosted GeoJSON data within QGIS. I have read a number of posts, such as add GeoJSON content in QGIS that shows very easily how to load GeoJSON into QGIS. I have had limited success with this process, however, for the most part I have not been able to consistently load hosted GeoJSON data into QGIS. I have tested with a few sources, including hosting my own data in Github with the same results (see the link below). I repeatedly get a not a valid data source error on load. My question, is there a standard for hosted GeoJSON data in order to load into QGIS, or are there settings that must be configured within QGIS?

https://github.com/ryangarnett/geodata/blob/master/locations.geojson


texts: PostgreSQL query using GROUP BY question

texts:
I have a polygon layer table in PostgreSQL  and I want to create a simple query using GROUP BY .

Here my fields :

fid   field1
--    ----
1     text1
1     text2
1     text3
2     text1
3     text4
4     text10
8     text1
1     text15
9     text20
9     text20

I want to use a PostGIS spatial query to group by my layer with fid field
but i need to keep all unique values from field1

Here an example where i want for export after query :

fid   field1
--    ----
1     text1,text2,text3
2     text1
3     text4
4     text10
8     text1
1     text15
9     text20

i try this Query :

SELECT fid,field1 FROM public.mytable GROUP BY fid

but I don't know how to keep all unique values from field1.

Any idea ?


texts: How to get monthly averages from Earth Engine in the Python API

texts:
I'm stumped.  It seems no one has done this in the python API.  I want to convert the Javascript below to Python.

 var sst = ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI')
 .select('sst')
 .filterDate(ee.Date('2013-01-01'), ee.Date('2017-12-31'))

var months = ee.List.sequence(1, 12);
var years = ee.List.sequence(2013, 2017);

var byMonthYear = ee.ImageCollection.fromImages(
  years.map(function(y) {
    return months.map(function (m) {
     return sst
      .filter(ee.Filter.calendarRange(y, y, 'year'))
      .filter(ee.Filter.calendarRange(m, m, 'month'))
      .mean()
      .set('month', m).set('year', y);
    });
 }).flatten());
print(byMonthYear)


texts: Attribute field based off intersecting polygon

texts:
I need to add a polygon for each bookmark (which works fine) and then attribute a field based off which Zip Codes the bookmark polygon intersects.

I don't know if I am on the right track. The code completes so there is no error to show you but the field for each polygon is empty.

cur = arcpy.da.InsertCursor(outFC, ["SHAPE@", "ZipCode"])
array = arcpy.Array()
for bkmk in arcpy.mapping.ListBookmarks(mxd):
    array.add(arcpy.Point(bkmk.extent.XMin, bkmk.extent.YMin))
    array.add(arcpy.Point(bkmk.extent.XMin, bkmk.extent.YMax))
    array.add(arcpy.Point(bkmk.extent.XMax, bkmk.extent.YMax))
    array.add(arcpy.Point(bkmk.extent.XMax, bkmk.extent.YMin))
    array.add(arcpy.Point(bkmk.extent.XMin, bkmk.extent.YMin))
    polygon = arcpy.Polygon(array)

    # Calculate Zip Code
    arcpy.MakeFeatureLayer_management(zipCodesSHP, "zipCodesLYR")
    arcpy.SelectLayerByLocation_management("zipCodesLYR", "INTERSECT", polygon, "", "NEW_SELECTION")
    zipCodes = [row[0] for row in arcpy.da.SearchCursor ("zipCodesLYR", "Zip_Code")]
    zipCode = ", ".join("\"" + str(zipCodes) + "\"")

    # Insert polygon and information
    cur.insertRow([polygon, zipCode])
    array.removeAll()


texts: How do I find out the resolution in (centimeters per pixel) of a TIF?

texts:
I tried using GDAL and got:

# gdalinfo nicholas_map_2.tif
Driver: GTiff/GeoTIFF
Files: nicholas_map_2.tif
       nicholas_map_2.tfw
Size is 28799, 27875
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-122.325695999999994,37.876846000000000)
Pixel Size = (0.000000441883000,-0.000000350276000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-122.3256960,  37.8768460) (122d19'32.51"W, 37d52'36.65"N)
Lower Left  (-122.3256960,  37.8670821) (122d19'32.51"W, 37d52' 1.50"N)
Upper Right (-122.3129702,  37.8768460) (122d18'46.69"W, 37d52'36.65"N)
Lower Right (-122.3129702,  37.8670821) (122d18'46.69"W, 37d52' 1.50"N)
Center      (-122.3193331,  37.8719640) (122d19' 9.60"W, 37d52'19.07"N)
Band 1 Block=28799x32 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=28799x32 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=28799x32 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=28799x32 Type=Byte, ColorInterp=Alpha

I assume pixel size is what I'm looking for? But is that in degrees?


texts: How do I toggle a legend of a cartodb layer in a leaflet map?

texts:
Per my previous question - here - how do I get the legend to toggle (I have it turned off now) in conjunction with toggling the cartodb layer via the leaflet layer control? Maybe adding it as a variable some way?

  var cdb_url = 'http://ovrdc.cartodb.com/api/v2/viz/52a89830-a8c1-11e4-b072-0e0c41326911/viz.json';
  cartodb.createLayer(map, cdb_url, {legends:false})
    .addTo(map)
    .on('done', function(crashes) {
      crashes.setZIndex(2);
      var overlayMaps = {
        "Vehicle Crashes": crashes
      };
      L.control.layers(baseMaps, overlayMaps, {collapsed: false}).addTo(map);
    })
    .on('error', function(err) {
      alert("some error occurred: " + err);
    });


texts: QGIS atlas generation to make company expansion slideshow/video

texts:
I need to create ~98 maps each showing the locations of customers as a company expands over a time series.

These will be made into a slideshow/video showcase so all maps need to have the exact same extent and scale regardless of customer expansion within (and beyond) the extent of the composer.

Each week's customers are currently in 98 separate CSV files.  I could get a single CSV with each week's batch of customers assigned a unique week ID if that is better.

I am unfortunately clueless about building expressions in print composer (and python) so my attempts to automate this process have all given me the wrong results.  I am very happy to share sample files, screenshots, (even the modest fee I am being paid to do this..) in order not to have to do it manually.

Ideally two boxes showing a) the week commencing date and b) the total number of customers active that week should be overlaid and updated automatically.


texts: How to center the text of a vector point feature in OpenLayers

texts:
I have several point features displayed in OpenLayers that, when a background fill is applied and small to no padding is set, don't center their text values in the center of the block fill:

This same text centers a lot better when the padding values are much larger (10)

The documentation states that for point features align, placement and baseline should have the text in the center of the block, but I'm forcing it just in case (this still doesn't work)

textStyle.setFont(`${style.FontWeight} ${style.FontSize}px ${style.FontFamily}`);
textStyle.setFill(new Fill({ color: `rgba(${ style.Colour.R }, ${ style.Colour.G },${ style.Colour.B },${style.Colour.A / 255})` }));
textStyle.setTextAlign("center");
textStyle.setTextBaseline("middle");
textStyle.setPlacement("point");
textStyle.setPadding([2, 2, 2, 2]) //[10, 10, 10, 10] centers the text way better

The font string value I'm sending in evaluates to: 300 9px Roboto

Is there anything I should be doing to improve the centering of these text labels?


texts:
Is it possible to find a map of Italian mobile phones antennas (and possibly also European ones)?

I need to import it in QGIS so it would be great a shapefile or dxf file or something related to them.


texts: How to programmatically enable / disable arrows control in Google Street View?

texts:
I'm building a web application using Google Street View andh HTML / Javascript.

I'd like to enable / disable the arrows controls that are on the images.

I can't show examples to do it ...

Suggestions / examples?


texts: How do I ignore points in a shapefile which do not overlay a raster in python?

texts: I have a geo-referenced tiff file and a point shapefile. I want to ignore all the points in the shapefile which are outside raster's boundary.


texts: Creating a single point on a specific distance on a line in QGIS 3.10

texts:
I'm trying to find a plugin of some sorts to help me create a point along a line at a certain distance. This does not require a import from a CSV file, just a field where I can insert a certain distance in meters. And where I can select the line where it should measure from and select from which side of the line I want to start measuring.

I'm searching for quite a while now and I'm only finding a lot of plugins which nearly do the same but most of these plugins require CSV import, or Chainage, which I don't need nor want.

I'm working with specialized electric fault meters which gives me distances along a cable where there should be a problem of some sorts, and it would be awesome if I could simply input the distance and get an exact location along the cable I've drawn.

PS: I'm not very skilled with coding, I can do some basics but that's about it.


texts: Reduce QGIS project needed PostgreSQL connections

texts:
I have a QGIS project, based on data from a PostgreSQL database. The project should be used by a lot people from my company. In a performance test I noticed, that every person, who opens the project block 3-5 connections to the database.

Is there any possibility to reduce the number of connections, which are blocked from the project?

The parameters, which I use for the database connection in QGIS.

connections\selected=QGIS
connections\qgis\service=
connections\qgis\host=123.45.67.89
connections\qgis\port=4125
connections\qgis\database=QGIS
connections\qgis\username=QGIS
connections\qgis\password=qgis
connections\qgis\authcfg=
connections\qgis\publicOnly=false
connections\qgis\geometryColumnsOnly=false
connections\qgis\dontResolveType=false
connections\qgis\allowGeometrylessTables=false
connections\qgis\sslmode=SslDisable
connections\qgis\saveUsername=true
connections\qgis\savePassword=true
connections\qgis\estimatedMetadata=false
connections\qgis\projectsInDatabase=false


texts: What information is needed to determine the section of the earth's surface visible in a photograph?

texts:
What information (either provided by a user or taken from sensor readings) would be required to determine the geographic area encompassed by a photograph? It seems like the area visible from a photograph would be a rectangle projected onto the surface of the Earth, but I'm unsure what information you'd need to figure out exactly what area a picture encompasses.

It seems clear that you would at least need the location, height, and orientation of the camera (but maybe I'm wrong about that?). There's probably also some information needed about the focal length, but I'm rapidly getting out of my depth here (pun regretfully intended).

As an alternative, suppose that you could get the user to geotag n points in the picture -- how many points would need to be tagged to compute the geographic bounds of the image?

The use case is that I would like to be able to present various statistics to a user about the area depicted in a photograph; ideally the user wouldn't have to do too much (e.g., if this were done from within a smartphone maybe all necessary information would be available from the phone's sensors), but if that's not possible it's okay to ask the user for necessary information.

Apologies in advance if this is too vague or posted in the wrong area; I'm happy to move it or add clarifying details as necessary.


texts: Merging multiple small overlapping polygons into one

texts:
I have data like this.

I need to combine these overlapping polygons into a single polygon. I used merging of polygons, but it takes more processing time in QGIS.

Is there any other way to do this?


texts: Can OpenLayers request different language tiles from map providers (eg Bing or Google)?

texts:
If I am in the US and my browser locale is en-US, I'll get English language maps by default from Bing.  Same thing if I access Bing tiles with OpenLayers.

If I change my country/region to (for example) Japan through http://www.bing.com/account/worldwide, then I'll see the entire Bing interface (including image tiles, place names, map features, UI components) etc in Japanese.

Is there a way for OpenLayers to request and render these fully translated tiles in a Bing or Google layer, based on a specific locale?  I'm primarily interested in these two providers, but if this is possible with OSM or others, even better!

To be clear, I'm not asking about i18n for the OpenLayers UI compenents such as the Layer Switcher, I'm asking about the source tiles themselves.


texts: Openlayers3 leaflet osm problem of displaying data

texts:
I'm a newbie in GIS world and recently built my own private osm tile server infrastructure that will be used in closed environment only and it works as a charm.

I tested the server inserting Leaflet code

var map = L.map('map').setView([48.1333,11.5667], 13);

L.tileLayer('http://{s}.mytileserver/{z}/{x}/{y}.png', {
    attribution: '&copy; OpenStreetMap contributors'
}).addTo(map);

L.marker([48.1333,11.5667]).addTo(map)
    .bindPopup('A pretty CSS3 popup.
 Easily customizable.')
    .openPopup();

The marker  exactly points  to Munich and it works perfectly.

Next I tried to do almost the same by using Openlayers example:

var centerpos = [48.1333,11.5667]; // position in degrees
var newpos = ol.proj.transform(centerpos,'EPSG:4326','EPSG:900913');
var map = new ol.Map({
    view : new ol.View({
        projection : 'EPSG:900913',
        center : newpos,
        zoom : 8.5,
        minZoom : 2,
        maxZoom : 12.5
    }),
    layers: [
        new ol.layer.Tile({
               source: new ol.source.OSM()
            })
         ],
   target : 'map'
});

and the map created shows somewhere near Djibouti - Africa.

What have I done wrong with Openlayers3 example ?

Or  did I do something wrong with installation of mapnik/postgis/osm ?


texts: Creating stable user Interface in MapBasic?

texts:
I am trying to have a user enter values via an interface and continue to run the program but the following line (Select...) fails to work most of the time as part of the program (although it works fine when you entered in a MapBasic window).

Dim Lat As Float
Dim Lon As Float
Dim sIDENT As String
Dim sValue As String

Control StaticText
  Position 150,50
  Title "Enter ID ?"

Control EditText
  Position 150,70
  Width 50 Height 15
  ID 99
  Value " " ' default value
  Into sValue

Select Lon, Lat from TableA where sIDENT = sValue into Querry


texts: How to import H3 Spatial Index modules into NodeJS App

texts:
I am trying to load a module called 'geojson2h3' in my NodeJS app. After running npm install geojson2h3, I tried following the initial code instructions given in Github to test it out:
import geojson2h3 from 'geojson2h3';

const polygon = {
  type: 'Feature',
  geometry: {
    type: 'Polygon',
     coordinates: [
      [
        [-122.47485823276713, 37.85878356045377],
        [-122.47504834087829, 37.86196795698972],
        [-122.47845104316997, 37.86010614563313],
        [-122.47485823276713, 37.85878356045377]
      ]
    ]
  }
};

const hexagons = geojson2h3.featureToH3Set(polygon, 10);
// -> ['8a2830855047fff', '8a2830855077fff', '8a283085505ffff', '8a283085506ffff']

const feature = geojson2h3.h3SetToFeature(hexagons);
// -> {type: 'Feature', properties: {}, geometry: {type: 'Polygon', coordinates: [...]}}

This method throws a SyntaxError: Cannot use import statement outside a module error. I then tried using the NodeJS require method as such: const geojson2h3 = require('geojson2h3');, replacing the import statement. This throws an Uncaught Error: undefinedModule error. The package shows in my package.json file, but I don't see it in my node_modules folder. I am assuming that's related to the issue. What did I do wrong here?
*Note: I know that this is a GIS-related package but the issue is a generic NodeJS package issue, not a GIS-specific one.


texts: Categorizing polygons instead of using default single symbol in ArcGIS for Desktop?

texts:
Is there any way to change display polygon from single symbol (default setting) to categorize?

So far, when I added a feature polygon into the table of contain, it will be displayed as only one symbol, but now I would like to have it always displays as category automatically without manually changing.

The default setting:

What it should be


texts: Leaflet Temporal animation of marker clusters

texts:
I recently have applied the markercluster leaflet plug-in to my map and I was wondering if something like leaflet.timeline would be compatible with cluster visualization?

if not is there any other temporal animation plug in that may work with the marker cluster within leaflet?


texts: Cannot Import Feature Class into Personal Geodatabase in ArcGIS 10.1

texts:
My objective is really simple. I have a shapefile and would like to convert it to a Feature Class. My goal is to use HEM tools, which I think requires me to convert my shapefile of interest to a feature class.

I have a geodatabase target where I have been able to create new feature classes and feature datasets just fine. Importing a feature class is really problematic though. I get different errors depending upon the method used. Right clicking on the geodatabase gives me the "Cannot Create Output Error" implying that the output folder cannot be found....but I can see it, and the tool can see it as well or it would let me run. Right-clicking on the shapefile and exporting to a geodatabase gives me the generic 99999 error. I've tried changing path names to eliminate spacing, have tried writing to different locations to avoid permission issues (their shouldn't be any since I am an administrator). Nothing works. I'm just missing whatever it is that I need to change.

Using the Feature Class to Feature Class tool works if I write to a folder...but then my end product is another shapefile, not a feature class. I new that would happen but just did it to test out writing to my Desktop Folder.

Oddly, I've tried writing to geodatabases that I have written to before and I can't import feature classes to them either. Strange since they are populated with feature classes that I have imported in the past.

I did just upgrade from 10.0 to 10.1...not sure if that has anything to do with it.

Thanks for your help,

Vin D'Angelo


texts: Extract shadows from a DEM using GrassGIS's r.sunmask command

texts:
I'm trying to use GrassGIS to extract the shadow from a DEM (elevation data).
The DEM I'm using shows Geneva (plus a part of the Alps mountains), and was captured by the SRTM. The elevation data was downloaded from the USGS website, and was already projected to the WGS84 coordinate system:

The command I've used in GrassGIS to calculate the shadow is r.sunmask, for which the date, time, and timezone has to be provided:
r.sunmask elevation=n46_e006_1arc_v3@PERMANENT output=alps-sunmask.tif year=2014 month=9 day=23 hour=11 minute=0 timezone=1

The problem is that the image produced with r.sunmask is completely empty. Is something wrong with the command above?
P.S: I've checked the answer suggested here, but I'm unable to make r.sun work since I can't find which parameter correspond to the output.


texts: Using QGIS processing tools in standalone script?

texts:
This is a follow on from Importing QGIS Processing in stand-alone python script?

I do not understand yet how to call an algorithm from the processing toolbox in a stand-alone Python script.

According to the Qgis documentation, I should use processing.runalg() but this doesn't seem to exist in my installation. When I do dir(Processing) I get:

['__doc__', '__module__', 'actions', 'addAlgListListener', 'addProvider', 'algs', 'contextMenuActions', 'fireAlgsListHasChanged', 'getAlgorithm', 'getAlgorithmFromFullName', 'getObject', 'getProviderFromName', 'initialize', 'listeners', 'loadActions', 'loadAlgorithms', 'loadContextMenuActions', 'loadFromProviders', 'modeler', 'providers', 'removeAlgListListener', 'removeProvider', 'runAlgorithm', 'runandload', 'tr', 'updateAlgsList', 'updateProviders']

Is the Qgis doc outdated or do I have a wrong version of something on my machine?

try:
    Processing.runAlgorithm('gdalwarp -overwrite -s_srs EPSG:32642 -t_srs "+proj=tmerc +lat_0=0 +lon_0=74.51666666666667 +k=1 +x_0=3300000 +y_0=14743.5 +ellps=GRS80 +units=m +no_defs" /home/steph/Temp/python_tmp/input/dsm.tif /home/steph/Temp/python_tmp/output/20160813-1225-UAVanalyse/dsm_proj.tif')
except:
    print 'ERROR ' + str(sys.exc_info())

Gives me:

ERROR (<type 'exceptions.TypeError'>, TypeError('runAlgorithm() takes at least 2 arguments (1 given)',), <traceback object at 0x7fa47898dbd8>)

But I'm not sure what the second argument should be?

While

try:
    Processing.runalg('gdalwarp -overwrite -s_srs EPSG:32642 -t_srs "+proj=tmerc +lat_0=0 +lon_0=74.51666666666667 +k=1 +x_0=3300000 +y_0=14743.5 +ellps=GRS80 +units=m +no_defs" /home/steph/Temp/python_tmp/input/dsm.tif /home/steph/Temp/python_tmp/output/20160813-1225-UAVanalyse/dsm_proj.tif')
except:
    print 'ERROR ' + str(sys.exc_info())

Gives me:

ERROR (<type 'exceptions.AttributeError'>, AttributeError("class Processing has no attribute 'runalg'",), <traceback object at 0x7f183586abd8>)

Using Qgis 2.16.1 on Ubuntu Xenial 16.04 LTS


texts: Adding shapefile causes ArcMap Drawing Errors (General function failure)?

texts:
I am creating a map that is a hodgepodge of databases, stand alone shapefiles, and a couple image services. Everything seems to go as expected until I add a certain shapefile to the mix. Then all of the other shapefiles disappear. Once that layer is removed from the map, everything else is fine. Only shapefiles are affected. The shapefile in question (as well as all others) display properly in 9.3. And, I'm adding other data that displays from the same folder, so it's not that I'm unable to access it. I've searched for similar incidents, but have came up empty.
A similar situation happens when moving the extent (10.7.1). I checked the projection and the geographic coordinate system. Everything matches the data frame. There are two File GDBs. error appears to only happen in the one that has data copied from a different GDB. The projection and geographic coordinate system are also the same for the GDB I'm copying from.


texts: How to customize the appearance of OSM maps on Garmin GPS device?

texts:
Just downloaded my first OSM map onto my Garmin etrex 30x from http://garmin.openstreetmap.nl/.
I always liked the extremly fine level of detail of OSM, as compared to other maps I know it contains really about every dirt path there is.

The map I just downloaded is working fine, however there are a few details I'd like to change.

There are way to many points of interest (grocery shops, gas stations, hairdressers etc.) marked on the map and I almost can't read anything. I'd like to delete these points category-wise (say, remove hairdressers but keep gas stations) to only keep what's relevant to me.
I'd like to change colors and line types of streets (say, make all trails a fine black dotted line and all highways thick and purple).
I'd like to adjust which elements show at which zoom level (say, only show trails when the area shown on screen is smaller than 20x20 meter).

Is there any software (preferably free) I can use to make these changes? Is there any website where I can customize these things before downloading? Or is there any other source of finely detailed maps I can use (relevant region is Central Europe)?


texts: Installing OSGeo4W on Windows, changed python version, "Import Error: No Module named site"

texts:
I've really been trying to install Proj4 on my windows 10 laptop, and make it accessbile via my git bash command. According to Proj4's documentation, using OS4GeoW is the easiest way to install the dependency (https://proj4.org/install.html) on windows 10 computers so I ran the setup.exe file.
I found additional instructions on proj4, however, that gave a command prompt instruction for installing proj4 from Os4Geo. The instruction was C:\Users\biney>osgeo4w-setup.exe --root C:\osgeo4w --quiet-mode --advanced --site  http://download.osgeo.org/osgeo4w --packages proj (biney is my name). And since I ran this on my command prompt... I got a pretty big output. But everytime I try running "python" or "python3," I get an error message saying "Import Error: No Module named site.py"

I was looking at this question and it seemed like a potential resolution was to set the pythonpath and pythonhome environment variables OSgeo4W64 Windows 7 & ImportError: No Module Named Site

I've done that through this pc --> properties --> advanced settings --> environemntal variables, whatever graphical interface that is provided... But now when I run python commands, I get the following error message:

 File "C:\Users\biney\AppData\Local\Programs\Python\Python36-32\lib\site.py", line 177
file=sys.stderr)
 SyntaxError: invalid syntax

This is after adding the PYTHONHOME and PYTHONPATH environmental variables. The command python --version still comes back as python 2.7.4. For clarification, the original and correct version of python I've been using on my machine lives in this file location: C:\Users\biney\AppData\Local\Programs\Python\Python36-32

I figured out the version of python switched by running python --version. It returns as python 2.7.4, but I can't quite verify which directory holds the python 2.7.4 version. My guess would be in OSGeo4W, but I could be wrong. This would be the file location C:\OSGeo4W\apps\Python36 (because there's also a folder in the same directory called Python27, but it doesn't have a .exe file).

I'm really stuck. How could I figure out where exactly the python .exe file that's being ran is being run from? And what do I do about this new syntax error message?


texts: Deleting field/column of shapefile with OGR/GDAL Python?

texts:
Is it possible to delete a field in a shapefile using GDAL with python?

With gdal 2.1.3 (Ubuntu 16.04), I can only see a CreateField() method on the layer object. No similar functions to delete a field.

Searching around, I can only find a ticket discussing the possibility and related rfc35.

Is this function implemented and how to use it?


texts: Selecting a subset of a raster, based on aspect and slope, as input to calculate area

texts:
Background
I created some digital elevation models and subsequently made a 2 layer raster for each containing the calculated aspect and slope using raster::terrain. Now I want to calculate the area of each which satisfies these two criteria:

135 < Aspect degrees < 215
Slope degrees > 2

Here is a link to a simple example of my dataset
https://www.dropbox.com/sh/5mfp421lg9jwkkh/AACG-_pTh56ie6zYvf6RRmRHa?dl=0

What I tried
I load both layers simultaneously as a stack. I first tried to do a simpler selection just to get a hang of the notation (being a beginner in R) on a sample dataset (I know its outside my own parameters set out above, but I just want to get it right first)

InvestigatedArea <- stack(paste0(getwd(),"/2150913_AspectSlope1.grd"))
result1 <- InvestigatedArea$slope[InvestigatedArea$slope[] ==  "2.336640"]
result2 <- InvestigatedArea$aspect[InvestigatedArea$aspect[] ==  "83.33418"]

Both these result in NULL answers, but I know my sample data has these values in them. Visible below in slots 81 in both instances.

> InvestigatedArea$slope[]
  [1]        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
 [11]        NA        NA        NA        NA        NA        NA        NA        NA 15.014699 27.780954
 [21] 28.800718 21.376579 14.421870  9.720741  9.799111 16.171215 22.605927 25.304127 23.419540 18.468254
 [31] 14.223759 11.670895  9.239475        NA        NA 24.906059 32.080647 27.778091 17.200335 11.943377
 [41] 12.618942 17.589796 24.648991 27.494242 26.727062 22.572401 15.775166 10.440066  8.643677  6.733919
 [51]        NA        NA 30.189089 30.811302 19.766258  8.661880  9.401229 16.495119 24.468115 28.694721
 [61] 25.355656 20.456161 14.836404  9.168124  5.176007  4.732281  4.400262        NA        NA 26.651682
 [71] 24.032007 12.866996  5.050277 10.976402 20.835245 29.398857 29.420214 19.465899 10.608719  4.032556
 [81]  2.336640  2.177207  2.326176  3.879496        NA        NA        NA        NA        NA        NA
 [91]        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[101]        NA        NA

> InvestigatedArea$aspect[]
  [1]        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
 [11]        NA        NA        NA        NA        NA        NA        NA        NA 307.59409 300.68588
 [21] 298.03110 296.04544 290.26505 293.71802 315.41357 322.72733 321.72003 318.79593 318.05820 319.33633
 [31] 316.50836 313.69717 320.50980        NA        NA 298.20007 303.81354 310.67441 317.35971 305.58899
 [41] 300.54056 309.42044 316.02017 322.52969 325.32404 329.24329 337.43942 330.18988 310.18207 309.40103
 [51]        NA        NA 298.71249 295.38638 302.32114 321.94833 303.58517 296.89182 301.38943 307.88290
 [61] 318.19205 322.84161 329.89163 358.89279 357.97546 296.30090 277.03336        NA        NA 299.59601
 [71] 284.44580 273.19974 277.92230 291.86438 292.40649 294.21609 300.87869 316.21036 322.60995 311.60651
 [81]  83.33418  89.59325 248.77850 236.84702        NA        NA        NA        NA        NA        NA
 [91]        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[101]        NA        NA

To be honest, I'm not even sure I am pointing to the right position where the actual points are stored. If I look at the structure of the rasterstack in the RStudio Environment browser, I cant actually identify the place where the data values are stored. I do know there is data, because I can plot it:

So I cant even begin to calculate specific areas, because I can't even get past this first hurdle. I know it's probably just a stupid notation mistake. What must I change in my notation to get a list of positions that contain the data that satisfies my conditions?


texts: web based GIS geoprocessing vs desktop based GIS geoprocessing

texts:
difference between the WEB and DESKTOP GIS geoprocessing (buffer)

what is the main difference in the area of processing speed or retrival speed of output
take in case of contour generation from DEM , what should I prefer in this case desktop or web based gis (in terms of time taken to complete task ,processing speed etc)


texts:
I have 250 road shapefiles, I have to append all the 250 shapefile into one single file. I have used ArcGIS append tool for combine layers (if number of shapefiles are minimum), but these 250 shapefiles take time as well as confuse me while selecting in ArcGIS append window.

Is there any other method in which I can append these 250 shapefiles easily?


texts: Error 000989 Python Syntax error in Field Caclulator

texts:
Using ArcMap 10.4 here and having some trouble running some python code in the field calculator.

I am trying to populate one text field based on the value of another text field. Here is what I am entering into the field calculator:

When I attempt to run this I get the error:

I am using 2 spaces for the 'if' line and 4 spaces for the 'return' line. Maybe I'm missing something big, brand new to coding here so any help would be great.


texts: Using ArcPy to ungroup group layer in MXD?

texts:
I am looking for a way to ungroup a group layer using arcpy... (instead of having to open the mxd and right-click on the layer to choose "Ungroup"...)

mxd = arcpy.mapping.MapDocument(pathToMyMxd)
df = mxd.activeDataFrame

newLayer = arcpy.mapping.Layer(pathToMyLayerFile)
arcpy.mapping.AddLayer(df, newLayer, "BOTTOM")

mxd.save()

layers = arcpy.mapping.ListLayers(df)
for layer in layers:
    if layer.isGroupLayer:
        # layer.ungroup()???


texts:
I have a layer of point data that includes fields for an ID number and an image file path. I would like to display the ID number for each point by default and be able to display the image in a callout for some points only.
To help illustrate the goal I have set up map tips that display the image using HTML. In the screen capture below the result is visible for point 17, but this only displays on the computer screen when a user hovers the mouse cursor over a point for a few seconds.
The label for point 21 on the other hand is displayed using a callout. I would like this callout for point 21 to do display similar to the map tip for point 17 so that I can display the image for specific points and have them visible on the final printed map. Is this possible to achieve?

For reference the HTML for the map tip is as follows:

[% BotanicName %]

[% right("Photo1",12) %]


texts: Using American Fact Finder to download data for all block groups in state at same time?

texts:
I am having issues using the advanced search on the American Fact Finder to download census block data for a whole state at a time. The best option that it gives me, when trying any route, is to select geographies of all census block groups within a county within a state. This leaves me with the option then of selecting every county, but in some states (like Texas) that have more than 250 counties, I am out of luck with this method.

Does anyone have a solution for this?

To duplicate my error:

Visit http://factfinder2.census.gov/faces/nav/jsf/pages/index.xhtml
Then select Advanced Search,
Then select on the left "Geographies",
Then select "all geographic types" radial button,
Then in the dropdown, select Block Group or part,
Then select state in the next drop down, and you will start to see
the issue.

Also, if instead of using the list search and you use the tab that says "name" I run into the same issue.

Does anyone have any experience with American Fact Finder downloads?

I am not asking the same question here as I am looking to see if this is possible just by accessing AFF not their ftp.
Gathering data for large number of census block groups?


texts: Obtain best-fit bearing for sequence of points

texts:
I'd like to calculate the overall bearing (direction) of a sequence of points in PostGIS.

It's clear how to do this for two points, using ST_Azimuth. How can I do this for the best-fit line through a sequence of points?

There is an example here on how to do it in ArcGIS, but is there a clean PostGIS solution?


texts: ArcGIS Pro matching and dissolving neighbours

texts:
I have a FeatureClass derived from ArcGIS Desktop 10.6.1's Polygon Neighbors tool [Windows 10 - 64]. I have supplied a tiny sample (of >250k) that gives the flavour and an example of the objective.
NEW LINK:
https://drive.google.com/file/d/1L58l02AWW64lhrE-lfZfLZ3PqXPlLf9K/view?usp=sharing
The attribute table is sorted by area - so we start with the smallest polygons and work our way incrementally through progressively larger features (there will be a cut-off beyond which I'll not want to progress - but let's not worry about that at this stage). In this example, we stop at OBJECTID = 7.
The attribute table includes fields src_OBJECTID and nbr_OBJECTID, generated originally by the Polygon Neighbours tool. I have dropped all neighbours other than the single feature with which the feature src_OBJECTID shares the longest boundary (attribute = Maximum_LENGTH).
What I am wanting to do is to dissolve progressively larger src_OBJECTIDs into their nbr_OBJECTID.
So, for example, the smallest polygon (src_OBJECTID '388') has nbr_OBJECTID '668'.
668 is thereby identified as the (longest shared boundary) neighbour of 388. 668 is also a (larger) SOURCE polygon feature further down in the table, with src_OBJECTID 668.
So I need do dissolve feature src_OBJECTID 388 into that neighbour feature - src_OBJECTID 668.
I then increment to the second-smallest feature. This is feature src_OBJECTID 131, which has to dissolve into its largest neighbour - feature src_OBJECTID 459. You will note that the third feature (OBJECTID = 3) is src_OBJECTID 459, and this has to be dissolved in turn into feature src_OBJECTID 668. So this is an example of nested dissolves, which are very common in the full dataset.
The subsequent features all dissolve into src_OBJECTID 668.
The difficulty I'm having is that I do not see how to select the target for the dissolve on the basis of it having the same src_OBJECTID as the nbr_OBJECTID of the fragment being dissolved into it.
Or, as @Hornbydd put it much more clearly:

... how to dissolve all the nested polygons, so in your example 131
dissolves into 459 which dissolves into 668 as well as all those that
dissolve directly into 668 which ultimately dissolves into 638?

I'll create an image tomorrow.
I am stuck at the first hurdle (no ModelBuilder stuff to show). I would first of all like to establish whether this is actually possible.
Then, if it is possible, establish the sequence of ArcGIS Pro tools required.