texts: Creating planar slices of 3D object to quantify voids?

texts:
We are trying to quantify the void space in a 3D object, very similar to measuring the holes in cube of swiss cheese.  I need to measure cross sectional areas and study the shape of the voids and we will conduct the analysis on multiple objects.  The cube is a multipatch object currently stored in an Esri geodatabase.  We have ArcGIS, Polytrans, and MeshLab software and are hoping for a no cost solution.

By no cost, I do mean beyond the cost of software that we already own.  We do have 3D Analyst, as well as Maya and 3ds Max.

I have been working with the "cube" file and have determined that it is not a closed object. I am new to 3D processing so I may not be using the proper terminology.

Here is an image of the cube in question.


texts: ArcGIS 10.1 Database Connection vs ArcSDE

texts:
I think one of the best capabilities of ArcGIS desktop 10.1 and 10.2 is Database Connection. Users can connect to every Database servers without using Middleware(ArcSDE) . You can Distribute, manage users  and other capabilities.But why ESRI release ArcSDE 10.1 or 10.2 while we can use Database connecion? What are capabilities of ArcSDE VS Database Connection?

Thanks For advance


texts: How can I create a valid probability map of species distributions based on presence/absence data from multiple quadrats?

texts:
I am somewhat new to GIS and only beginning to learn spatial statistics/analysis. I am an ecology grad student doing a before and after study on invasive plants and their responses to management efforts. First, Here are a few maps to clarify my design and avoid confusion.

My experimental design consists of ~100 transects scattered throughout a park. Each of these transects has anywhere from 1 to 10 quadrats placed 5 meters apart along its length, and each quadrat has an area of 1 m^2. Within each quadrat, the presence/absence of 11 different invasive species was documented, first as baseline data in 2007, and again in 2014 to determine if management efforts in the park are controlling their spread. However, some areas of the park have transects that are placed very close together (25 m apart) while other areas have transects placed 50 m apart, and I am concerned about how to deal with the possible oversampling and spatial bias in these closely placed transects. I had no control over the original design and baseline data collection, and to me this design is just plain ugly.

Based on the presence/absence data I have collected, I would like to know each species' distribution for each year (22 maps), and be able to determine if there are any statistically significant hotspots of certain invasive species, or if certain species are correlated with one another. I had ArcGIS calculate the centroid for each quadrat so that I have both polygon and point features to represent the quadrat, if need be.  So far, I have created several Kernel Density maps using the centroid as my point data and the column with binary data (presence=1, absence=0) as my population field, and leaving the search radius blank, since I am unsure of what to do about the variable distance between transects. I'm almost certain that Kernel Density is not what I should be using to generate this map based on what I've read so far on here.

Here is an example of one Kernel Density for Euonymus fortuneii in 2007:

Qualitatively, it does appear that these maps match up with what I casually observe in the field, but I need raster output files that do not violate any key assumptions and adhere to proper modeling methodology. I need an output raster that I can perform quantitative stats on and get information on significance. From what I have read, it seems like interpolation as some type of kriging is what I need to do, but I am not familiar with these techniques. However, I do have the license to run all the necessary extensions in ArcGIS and can (hopefully) learn them if I have to.

I would also like to create change detection maps for each species, showing areas that their presence has decreased and areas that it has increased, and whether or not these are significant. I have used the raster calculator tool to subtract the 2007 Kernel density rasters from the 2014 rasters to produce a map that shows areas of decrease in blue, and areas of increase in red.

I can tell based on spreadsheet data alone that Alliaria petiolata has significantly increased in presence in 2014; it was found in 76/464 quadrats in 2007, and 176/464 quadrats in 2014. I would like to be able to do this in ArcGIS with raster files though, and I know there has to be a way.


texts: Remote SQL Server, Crystal Reports and Geospatial Join

texts:
I have:

A remote SQL server (2008 R2) that I have access to (But NOT full control of) which contains my source data including OSGB 1936 Geometry [ESPG: 27700]
In house SQL Server 2008 R2
A Crystal Report Automator in house that allows me to send multiple variants of a report in PDF format to multiple recipients (500+)

Crystal XI

limited comprehension of CLR Package installations!

Because of an anomalous data rationale during the build phase I have a need to connect two tables of spatial data in order to get a full realisation of the required report.

I can't get this done on the on the remote server for a variety of logistical reasons.

I had attempted to build Views (of the simplest Select * from mytable variety) that were giving me a degree of success in SQL when I invoked the OPENQUERY option, which had functioned but now are returning me:

The OLE DB provider "SQLNCLI10" for linked server "xx.xxx.xx.xxx\MSSQL2008R2" supplied inconsistent metadata for a column. The column "Geom" (compile-time ordinal 4) of object "SELECT ColimnID, Column2, Column3, Geom, Column4, Column5, Column6, Column7, Column8, Column9, Column10
FROM remotedata.dbo.mytable
 " was reported to have a "Incomplete schema-error logic." of 0 at compile time and 0 at run time.

In the absence of being able to perform the joins within SQL directly because of this I'm instead trying to do some of the geo-manipulation in the Crystal context.
This however is presenting a secondary set of challenges;
what could ordinarily be expressed as

mytable.geom.STX

now needs to be

geometry::STGeomFromText(CONVERT(nvarchar(50),"mytable"."Geom"), 0).STX

or at least in that construct I can get it to work when entered as a SQL expression field within Crystal.

Unfortunately the really thorny Join that I need to perform

mytable as a INNER JOIN
mybasereference AS b ON b.geom.STContains(a.Geom)=1

and I have absolutely come to a dead end as to how (or even IF) this is theoretically possible within Crystal let alone how to make it read correctly for Crystal to process.

I'm agnostic as to where I invoke this join within the entire processing as it only provides me with a reference- but a necessary one

Is there any way that I can overcome the incomplete Schema-Error logic error [and still perform the join that I need]? (The remote server controller has informed me that they have made no changes on their side and that it must be an issue with my server)
Is there any way to perform a spatial join like this in Crystal?


texts:
I'm trying to get the user I used to login to Portal in ArcGIS Pro, but I can't find how. Is it possible?

I manage to fetch the Portal information using

arcpy.GetPortalDescription()

, but not the user used to login in.


texts: I'm new to PostGIS and pgRouting and I'm using the pgr_dijkstra functionality and I need to convert the values it returns in the format lat/long, is that possible? If so, how?


texts: Opening geoprocessing tools causes ArcGIS to hang

texts:
I am getting the same problem as http://forums.arcgis.com/threads/19468-Opening-tools-causes-ArcMap-to-freeze-crash

That is, when I open any geoprocessing tool in ArcCatalog (might happen in ArcMap as well, though it hasn't to me yet), sometimes the application hangs. The timing of the hang occurs after the tool dialog opens but before it finishes drawing. I don't know exactly what conditions cause this to happen, but it feels random, maybe in 1 in 4 attempts it hangs.

For me, terminating the hung process and trying again sometimes works, sometimes it doesn't. But it's a huge headache.

Has anyone come across the same problem and found any resolution? I have already tried reinstalling ArcGIS. I think this started happening some time after installing SP3 for ArcGIS 10.


texts: Merge neighboring image tiles into a single image

texts:
I have a set of satellite image tiles from Bing Maps (info here). The tiles' names specify tile coordinates within the tile system and zoom, e.g. 156608-142775-18.jpg.
An example set of tiles:
['156608-142775-18.jpg', '156608-142776-18.jpg', '156607-142777-18.jpg', '156608-142772-18.jpg', '156607-142776-18.jpg', '156608-142777-18.jpg']

Here are the tiles visualized:

What I need to do is to merge all neighboring images into one image. So, in my example, I wish to merge the lower five tiles into one and output two images - the big merged image and the upper single tile which has no neighbors.
I am not sure how to achieve this. It seems to me there should be a simple solution to this that I cannot figure out.
I will use the same code for merging a large number of tiles so I want to avoid solutions that use a lot of memory or take long time. I work in Python. I don't need the actual code, just an idea how to do this.


texts: How to read this data file and convert it to raster format for analysis?

texts:
I am trying to introduce the 5 degree dataset from here http://www.cru.uea.ac.uk/cru/data/precip/ (it is the first link). And once I introduce it to GIS to use it as a raster input to calculate zonal statistics using as feature zone the world borders (http://thematicmapping.org/downloads/world_borders.php).

My first question is how to introduce the dat data into GIS.
An second, once I introduce it can I use it as a regular raster file and do zonal statistics?

Can anyone help? I use ARCGIS 10.1

Thanks!


texts: Shapefile failing to import using shp2pgsql?

texts:
I'm trying to import a shapefile (boundary.shp) using shp2pgsql. I'm using the following query

shp2pgsql -I -s 4326 boundary ch01.test_boundary > psql -h 127.0.0.1 -d spatial_database -U postgres

After running the above import statement I get the following message:

Shapefile type: Polygon
Postgis type: MULTIPOLYGON[2]

When I query for the table ch01.test_boundary it wasn't created, so it's failing.

I can successfully import the same shapefiles (get the same type messages, but the data is imported) using pgAdmin so I suspect it may have something to do with permissions.

I did notice i am not prompted for my password when running the above statement.

How can i resolve this?


texts: How to visualise change matrices over time

texts:
i've hit a bit of a tricky visualisation issue;

I am quantifying some land use and land use change (LULUC) over a number of years. Typically, for data going from year to the next, a land change matrix is pretty handy;

In the above image you can see what percentage of the total land present in 2005 (the left column) becomes/stays as what in 2006 (the row along the top). Eg category A staying as category A is 50.2% of all land in the dataset, category B changing to category F is 2.4% of the land etc.

it is important to realise that this represents quantities of change from 05 to 06, not absolute quantities of that land category in either 2005 or 2006.

However i have a dozen change matrices, all important, but i'm slightly loathe to show them all in my results.

I thought this multivariate time series idea for R was pretty nice, but it will only show absolute quantity from one category to another.
A lattice/trellis plot where by the land categories are the 'outer' axes forming the structure of the trellis while 49 bar charts showing % change over each year step would be displayed, but this sounds horrid.
A 3D plot where by the matrix above is in a 3D plot like a set of tiles, the colour of the square is used as the quantity indicator (think heatmaps and geom_tile in R), and then the data in the time series is stacked on top of it, almost like a stacked bar chart but each bar would consist of 'tiles' reaching the same height but colour to denote quantity at each level. Sounds complicated.

Are there any good ideas out there for this issue? perhaps multiple matrices are the clearest way.

(p.s. R is preferable but anything that works is fine too :)  )


texts:
I'm new to pyqgis and trying to bring up from the previous post with a little progress from the last (implemented debugging tool pdb) but now trying to export into a .tif file so that I can check on the vertices from the csv geometry. Unfortunately I get stuck with an blank render .tif image. The render.tif filesize is showing 1.9mb, when I download and view locally, the image only has the background with white colour and no vertices.

At the same time, I checked my csv file and also looping through the features of the layer after export the csv files. They are all good...

How can I have all the features from the vector layer to be shown on the .tif image file?
Anyone please help? :/

PS: I'm on ubuntu server and using PyQGIS.

Tried with render.fullExtent() and render.extent(), both return blank image.tif

Inserted sample of elevp.csv

My code:-

#!/usr/bin/python
import sys
import qgis
import PyQt4

from qgis.core import *
from qgis.utils import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtXml import *
import pdb

qgis_prefix="/usr"
QgsApplication.setPrefixPath(qgis_prefix, True)
QgsApplication.initQgis()
app = QgsApplication([], False)

uri = "/home/ubuntu/pyqgis/elevp.csv?crs=%s&delimiter=%s&xField=%s&yField=%s&elevField=%s" % ("EPSG:4326",";","x","y","elev")
vlayer = QgsVectorLayer(uri, "test", "delimitedtext")

QgsMapLayerRegistry.instance().addMapLayer(vlayer)

# create image
img = QImage(QSize(800, 600), QImage.Format_ARGB32_Premultiplied)

# set image's background color
color = QColor(255, 255, 255)
img.fill(color.rgb())

# create painter
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)

render = QgsMapRenderer()

# set layer set
lst = [vlayer.id()]
print(lst)
render.setLayerSet(lst)

# set extent
rect = QgsRectangle(render.fullExtent())
rect.scale(1.1)
render.setExtent(rect)

# set output size
render.setOutputSize(img.size(), img.logicalDpiX())

# do the rendering
render.render(p)

p.end()
# save image
img.save("/home/ubuntu/pyqgis/render.tif","tif")

Sample elevp.csv:-

x;y;elev
120;960;13
360;40;52
640;840;3


texts: Add Date to raster layers in raster stack in R

texts:
I have a raster stack and I want to add to it a Date field which will include the year-month-day for every rasterlayer. I have my dates in a one column data.frame.
I tried to do this with addLayer but it didn't work.

Any ideas?

Example:

dt<-as.data.frame(as.Date(c('2010-11-1','2008-3-25','2007-3-14','2007-3-18','2007-3-20')))

ras = raster(matrix(1:30,3,4))
my_s = stack(ras,ras,ras,ras,ras)
my_s <- addLayer(my_s, dt)


texts: How to create circle that contains a certain number of points?

texts:
In ArcGIS 10.2 (with spatial analyst ext). I have two point layers.

One represents schools and one represents students. For each school I have a number of students living near the school I need to "grab" with a circle. The circles (a polygon layer) are something I will need later for another process.

I could have sworn ESRI has a tool to do this but I can't find it or remember key words to aid my search.

Example of what I am looking for: School(point) has an attribute of 100. I also have a point layer with thousands of students. I want to create a circle polygon (with school as the center) that contains 100 students.


texts: Arcpy UpdateCursor not Updating, but no Error Reported

texts:
Below is the code I'm using in an SDE GDB, on ArcMap 10.4.1:

import arcpy

# Variables
inputAsbuilt = arcpy.GetParameterAsText(0)
transmedia = inputAsbuilt + "\\Transmedia"
structure = inputAsbuilt + "\\Structure"
spliceClosure = inputAsbuilt + "\\SpliceClosure"
span = inputAsbuilt + "\\Span"
equipment = inputAsbuilt + "\\Equipment"
buildingfootPrint = inputAsbuilt + "\\BuildingFootPrint"
attachment = inputAsbuilt + "\\Attachment"

spliceClosure1 = '\\\\Bvweb104draco\\e\\03_MapServices\\SDE_CONNECTIONS\\bvtampadev.sde\\bvtampadev.sde.Telecom\\bvtampadev.sde.SPLICECLOSURE'
buildingfootPrint1 = '\\\\Bvweb104draco\\e\\03_MapServices\\SDE_CONNECTIONS\\bvtampadev.sde\\bvtampadev.sde.Telecom\\bvtampadev.sde.BUILDINGOUTLINE'

fcList = [transmedia, structure, spliceClosure, span, equipment, buildingfootPrint, attachment]
sdefcList = [spliceClosure1, buildingfootPrint1, fiberEquipment1, conduit1, structure1, fiberCable1]

# Print out FQNID's of feature classes inside As-Built GDB
fqnList = []
for fc in fcList:
    with arcpy.da.SearchCursor(fc, "FQN_ID") as cursor:
        for row in cursor:
             fqnList.append(row)

for fqn in fqnList:
arcpy.AddMessage(fqn)

# Update "bv_asbuilt_gdb_uploaded" field in SDE feature classes
for fc in sdefcList:
    with arcpy.da.UpdateCursor(fc, ["bv_fqn_id", "bv_asbuilt_gdb_uploaded"])      as cursor:
    for row in cursor:
        if row[0] in fqnList:
            row[1] = "Y"
        else:
            row[1] = "N"
            cursor.updateRow(row)`

What I want the code to do:

Make a list of unique values from a field in the input feature class
Take that list of unique values, check it against a field in another feature class to see if they're ==, if not write N in the "FQN_ID" field, if yes wrote Y in the "FQN_ID" field

The code is throwing no errors and appends a string "N" where I want it to, but it's not checking the row[0] against "fqnList" correctly, so it populates the entire field with "N" even if there the list matches a value in the feature class field and should populate as "Y"


texts:
I am planning to go on multiple days hike and used Wikiloc to download GPX tracks. We will use a GPS for our navigation but I always want a printed map as a backup.

Maps for that area are notoriously bad but OpenStreetMap contains most of our trails and Points of Interest.

I loaded the track in QGIS with the OpenStreetMap and this is a good start. The only thing that is missing are contour lines (elevation lines) and potentially hill-shading.

What is the easiest way I could add contour lines and hill-shading to my map (for a printout)? Ideally, I would add a WMS or VectorTile layer and not go through the process of downloading a DEM and creating everything myself.

UPDATE:

It is fairly easy to add a mapbox WMTS layer to QGIS. I opted for the Outdoors mapbox style which has both contour lines, hillshading and OSM trails. However this introduced another problem. When printing / exporting the scale of the WMTS layer changes and becomes unreadable.

However when exporting using the print composer this is the result:

similar issue
another


texts: Using ModelBuilder to update feature line length

texts:
I maintain a network of lines (geometric network) in an enterprise geodatabase and would like to automate the updating of a field after I complete my editing. The goal is to run the script nightly on a server to update the field.
The field I'd like to update represents a custom length value (miles) that is different then the shape length value.
Is it possible to create a model using ModelBuilder that will allow the calculation and updating of this field?
What tools would you use in ModelBuilder to calculate geometry and update the rows?
It does not matter if all rows are recalculated nightly even if they are not edited.
This is my first attempt in trying to build a model that addresses this type of update.


texts: Creating buffer around point in LeafletJS?

texts:
I am using LeafletJS. I want to create buffer around point.

After Selecting a point buffer radius will be specifed by user in feet.
After searching of a while I found this link of JSTS.

Does JSTS buffer works with LeafletJS? I am not able to find any sample for this. Any other solutions to achieve same will be appreciated.

After looking into the code of JSTS and searching for a while, I came to a conclusion, there is no conversion between JSTS geometry and Leaflet JS geometry. So it is not possible to use it without implementing conversion.


texts: How to remove banding in Sentinel-2 imagery over water?

texts:
In several L2A Sentinel-2 images, I've seen an angular banding effect over water bodies. Another user thought that these seam effects were caused by aerosol scattering affecting each of the S2 MSI detectors. In any case it would obviously be nice to remove this type of artefact. My ideas right now are

Rebuild the L2A product from L1B using different parameters in Sen2Cor. Maybe the wrong default parameters were used in the atmospheric correction?
Use a SNAP tool or plug-in that I don't know about yet.
Write some custom software to bias correct each detector.

There are some suggestions in this paper, I just don't know where to start or what the current recommended solution would be.


texts:
I am new to QGIS and just downloaded QGIS 2.4 Chugiak.

I am practicing QGIS by using training manual, as mentioned in above subject.

I have questions in regard with the sentences in page 7:

  You should then extract and save the following fields from the lines
  and points layers to their corresponding directories:

  lines: “highway != ‘NULL”’ to roads, and “waterway != ‘NULL”’ to rivers
  points: “place != ‘NULL”’ to places

How do I extract and save?

Can someone help me?


texts: Mean of polygons area within another polygon layer

texts:
I have the following layers :

Sectors (polygons limited by the black lines)
Houses (polygons in light blue)

I have created the area attribute in the Houses table. I have also created an attribute that represents the centroid of the geometry.

I would like to calculate the mean area of Houses located in each sector and add it as an attribute of the Sectors table.
A house is considered "in the sector" if its centroid is in the sector (to avoid problems with houses over 2 different sectors).

I don't know how to do this "group by" operation and haven't found answers searching this stack.


texts: Setting transparency for range of colors in raster layers in QGIS

texts:
The starting point: Loading a raster to QGIS, you can set individual colors to transparent in the Panel Layer styling / Tab Transparency by activating the icon Add values from display and click on the map canvas. The color of the selected raster pixel will be set to transparent (see screenshot below).
The problem and what I tried: As can be seen in the screenshot, often (especially using georeferenced jpg-images) what for the eye is one color, technically consists of several different colors with very similar values for red/green/blue: different shades of the same color. Of course, you could repeatedly use the Add values from display and click on the map to add further colors to the transparency setting.
The question: However, this can be annoying. Is there a way to select a range of colors to be set to transparent? I'm looking for something like tools known from image processing software like Photoshop or Affinity Photo, e.g. the Flood Select Tool in Affinity Photo: Tolerance—sets the range of pixels affected when a pixel is clicked.
Screenshot: Only part of the see is set to transparent (RGB value: 221/242/247), showing the OpenStreetMap basemap in the background. Neighboring pixels with - let's say - RGB 225/240/243 are still opaque:


texts:
I used extract() function to get remote sensing image cell value.The same code to process tif file produce error message. Code as follows：

r1 <- raster("2000_2605_1.1_km_monthly_NDVI.tif")
polygon_mus <- extent(c(107, 111, 37, 40))
cell <- extract(r1, polygon_mus)

This *.tif is a output file from MODI hdf file with MRT. The process result is as follow:

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function extent?for signature ?NULL"?

Where is the error? I seem to be missing it.


texts: Is it possible to mix italic and normal font in a label?

texts:
I have concatenate the values of two field like

"BYKV_BTYP"  || '/' ||  "FFH"

Is it possible to draw the value of the column  "FFH" in italic if it is not empty and the value of "BYKV_BTYP"  normal? Such as

G212-6510


texts:
I have remote sensing data which I have got from landsat 7. The area which I am working on is combination of farms and natural area (Bare earth; without any vegetation). I want to estimate the Evatransporation for whole area to put in my hydrological model.

Can I use NDVI index for this purpose?

As I know NDVI helps us to find evatransporation in the vegetated area.

How can I estimate the evaporation from bare earth regions?


texts:
I am using the GRASS Python Scripting Library on Ubuntu 14.04 with the packages available from the Universe repository. I have the following block of code:

        try:
            grass.mapcalc("$weightLine = $seg_rast * $cost",
                        weightLine = "weightLine",
                        seg_rast = "seg_rast",
                        cost = "cost")
        except:
            self.logger.debug("Exception running mapcalc: " +
                              str(sys.exc_info()))

        kwargs = {}
        kwargs['stdout'] = grass.PIPE
        kwargs['stderr'] = grass.PIPE
        ps = grass.start_command("r.sum", rast="weightLine", **kwargs)
        mesg = ps.communicate()
        self.logger.info("r.sum: %s", str(mesg))

Which is logging the following:

[2014-06-03 09:57:35,532] DEBUG: Exception running mapcalc: (<type 'exceptions.SystemExit'>, SystemExit(1,), <traceback object at 0x7ffcbab1de18>)
[2014-06-03 09:57:35,538] INFO: r.sum: ('', 'ERROR: Raster map <weightLine> not found\n')

The employment of sys.exit() is not very helpful, even though it can be worked around. Without the except: block it exits my script and any other processes that may have invoked it. Note also that except Exception, e: does not catch exceptions.SystemExit.

But worst of all, the mapcalc instruction is not succeeding. It is either failing - in which case an error message is lacking - or it is simply not being executed.

How can I know the cause of this?

Update I: I have tried to run grass.mapcalc in different ways:

grass.mapcalc("weightLine = seg_rast * cost")

And:

weightLine = "weightLine"
seg_rast = "seg_rast"
cost = "cost"
grass.mapcalc("%s = %s * %s" % (weightLine, seg_rast, cost))

But the result is always the same: exceptions.SystemExit and the command is not executed.


texts: Selecting the name of the raster layer for each point in a shapefile (QGIS)

texts:
I have multiple raster files (i.e. elevation in the United States) and over 300 points in a separate point shapefile. Each point lies in only one raster file. I'd like to be able to find out the name of the raster file that each point lies in.

Ultimately I'm trying to do a batch run over all the points, but to be able to do that I need an input raster file for each point hence why I need a way to associate  the name of a raster file to each point

Does anyone know an easy way to do this in QGIS?


texts: We are moving to France soon, and are currently investigating a lot of areas.  Does anyone know of good online/offline resources for French mapping?  I am very interested in mapping, weather, environmental, geological, hydrological data, as we are preparing to create a zero emissions home out there and want to get the best area.


texts: Changing attribute length of existing field afterwards?

texts:
I already have about 100 values in a field in QGIS. I set it to 50 characters because I was thinking that would be enough. But that was wrong. I need more space in the column. But there is no way to change it. In the usual editor you only can remove or add a new field.

And in the "advanced" table manager you can change the position and rename it. But nowhere where I edit the length of a field.

The only thing I could do is creating a new field with more space. Then import the data somehow. But does that work somehow? Or is there a way to edit the text length afterwards. I'm also interested in hacks.

And yes I really need that. There is no way to get the string any shorter. I think 100 chars are fine.

Now it came to mind the following. Could maybe someone tell if that is working before I try that. So First of all I would create a new field with 100 chars column. Than select all old features and CTRL-X them and CTRL-V into the new shape.

Will that work out?


texts: Get image pixel centroids in Google Earth Engine

texts:
I would like to create a feature collection consisting of the centroid points of the pixels in an image. The question of how to do this has been asked before on Stack Exchange, but not actually answered.
For example, suppose we have a single image such as the ESA 2009 global land cover. I crop this to the boder of the Republic of Congo:
var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');
var congo = countries.filter(ee.Filter.eq('country_na', 'Rep of the Congo'));

var dataset = ee.Image('ESA/GLOBCOVER_L4_200901_200912_V2_3').clip(congo);

With this data, how do I obtain a collection of points at the centroid of each pixel in the image? Converting the image to a feature collection does not seem to be the right approach, as this assigns polygons based on properties (not pixels).


texts: Finding and copying text with Field Calculator in ArcGIS for Desktop?

texts:
I'm trying to find every instances of an "Avenue #" in the Name field and have Field Calculator copy that over the single letter written in the Subname field so that both fields look exactly the same, but neither python nor visual basic are a strong suit of mine.

I need ONLY the avenue (letter) texts copied. Is there a way to do this with field calculator?

EDIT: Adding this for clarity...

I am needing only the "AVENUE #" text within the NAME field to be copied over to the SUBNAME "#" fields. The results would be an example such as this:

SUBNAME "C" = NAME "Avenue C". But SUBNAME "Avenue Belle Chase" would remain untouched.


texts:
This is in ArcGIS Pro.
I will be performing a raster calculation to determine how much a given raster might change over time by multiplying the current raster by the percent projected change and then adding back to the current raster to give the 'new' raster
The problem is that the values for the raster need to be constrained to -17 (the worst case) to +17 and this approach could make me fall outside of those values.  So, I know I need something that works sort of like an if then statement - if >17 then =17 and if < - 17 = -17  This is the step that is tripping me up.  I have never been able to really come to grips with the CON statement and I think this may be what is needed here so could you provide any examples?


texts: Reording fields in Identify Results panel

texts:
What can be done to modify the fields order in the Identify Result panel, idealy without re-creating a new layer?

Modifying the table field order or the UI field order, as described in this answer, does not affect the  Identify Result panel.

Having a layer with fields a, b and c, I would like to show c,b,a (or any other custom order). I am currently using QGIS 3.4.15 LTR.


texts: Selecting the longest connected lines from PostGIS ST_ApproximateMedialAxis

texts:
I used PostGIS ST_ApproximateMedialAxis to transform polygons representing mountain ranges into a series of lines.  Now, I'd like to select only the "longest" line segment for each group, so that I can dissolve the lines and use it for "curved" labeling.  Can I extend the query in PostGIS to end up with the single series of arcs?

The end result should look something like the lines highlighted in yellow in the picture below.


texts: Cursor to update feature class by location?

texts:
I am trying to update an existing field in a polygon fc with an existing field in a point fc. There is only ever one point fc located within each polygon fc. The field names in the polygon fc can not change as it needs to have a matching schema for a later append.

Lotnumber is the point fc while fcUpdate is the polygon fc in the below code. I was trying to create a search cursor that would iterate through each point and store the value in the TEXTSTRING field. the code should then select the polygon that contains the current point, insert an update cursor, and update the Housenumber field with the value obtained from the point.

# Enter the point feature class containing lot numbers
Lotnumber = arcpy.GetParameterAsText(4)

# Make a layer from previously created new feature class
arcpy.MakeFeatureLayer_management(fcUpdate, "Parcels_Temp")

# open a search cursor to iterate through each point
with arcpy.SearchCursor(Lotnumber, "TEXTSTRING") as cursor1:
    # Iterate through each point in the point fc
    for row in cursor1:
        lotnum = row.getValue("TEXTSTRING")
        # Select the parcel that contains the point
        parcel = arcpy.SelectLayerByLocation_management(fcUpdate, "COMPLETELY_CONTAINS", row)
        cursor2 = arcpy.da.UpdateCursor(parcel, 'HouseNumber')
        for row in cursor2:
            row[0] = lotnum

# delete the cursor
del cursor

My code is based on a question answered here: Using Select By Location to update field in feature class using ArcPy?

I'm very new to Python. I've tried running the script in ArcMap and get an error stating invalid SQL statements were used.


texts: Why does ArcMap tell me that I must be signed in to ArcGIS Online when I use the Find tool? And how do I take this setting off?

texts:
I recently upgraded to ArcGIS 10.4.1 and overall I am liking it very well. There's just one annoying problem that has been occurring ever since the upgrade.

When I use the 'Find' tool after it has been closed out for awhile, ArcMap seems to hesitate and then pops up with the message: "You must be signed in to ArcGIS Online with an organization account to use the ArcGIS Online World Routing Service."

Here are the exact steps I follow:

open find features window

input account number
select layer to search and attribute to search
search and find result
right click result (this is when it seems to take longer than it used to)
select 'zoom to' in drop down menu
ArcMap zooms to property, blinks it and then instantly displays the message
generally I hit cancel and continue my work, but it has become a nuisance

I definitely have a Legit copy of ArcMap and I do have an organization account at ESRI, but I don't want to be told to log into it everytime I search a property. I am not running any ESRI layers nor searching any world files. I am searching localized data that we build in house.

Why is this happening, and how do I fix it to where I'm not bombarded by this message everytime I search?

Edit: Adding screenshots of the other two tabs in the find tool for reference.

It only occurs when I right click the items in the results field. If I just do a quick double click it will zoom to them without displaying to message. It's just difficult for me to do that, because the way I learned was to right click and select zoom to.

Also it only occurs like the first time I use the find tool after it's been closed for awhile. So it's not absolutely preventing me from getting work done, I just know with ArcMap that almost everything comes down to a hidden checkbox or something that I'm not finding.

Edit: Adding additional screens to show this occurring with multiple layers within the same map.

I cannot, however, seem to replicate this with some of my more recent maps that I've created since updating. As per Simon's comment I will create a fresh version of the same map and see if that is the solution.


texts: Snapping two topology layers together based on point

texts:
In QGIS 3.14.15 with GRASS 7.8.3 I have two layers both with correct topology that function how they should. When using the Network Analysis > Service Area tool, they work as they are supposed to. Trace and Shortest path etc.
However, I would like to connect them to each other.
What it currently looks like:

When running a trace on the RED it should also expand the network to the BLUE as they 'SHOULD' be connected at the point layer.
I know this is possible with manually editing the topology and simply dragging it to the point or to each other. This works well,  though this network is way to large to do manually and I need an automated tool.
I have used grass v.clean to clean up each layer individually, though when trying to use it to connect the two it changes the network too much as it uses all nodes when using the 'snap' tool in v.clean. The distances between the network layer(s) and the point layer range between .5m and 2.5m
What I want it too look like: (as mentioned it is possible manually)

Is there any 'black hole' or something that would just 'suck' all points to the point layer? Slight variations in the network are acceptable, such as when I adjust it manually. Just need to do this for a couple million points/lines.


texts: How can i use openModeller plug-in for QGIS?

texts:
As far as I know, niche modeling is possible in the old version of QGIS by using openmedeller plug-in.
I want to use 'openModeller' using the current version of QGIS(v 1.7.4).
Is there any solution to make it work?

I would appreciate your help.


texts: Feature counts and rule-based styling for overlapping points of diff. height with one height class per map

texts:
I have a set of 700 points within the footprint of a 16-storey building. The points have a building floor attribute (e.g. "Floor" = 12), and another attribute relevant to my presentation (call it "AttType" = 'A', 'B', 'C'...).
I want to generate a map for each floor which should do two things:

display only the points on that floor
show accurate point counts per AttType in the legend

I have managed to do 1 by creating 16 map composers and naming the map ID "Level 1", "Level 2", etc., then styling the point layer with the rule @map_id = concat('Level ',"Floor") so that only points of the appropriate floor show up.
However, because of this non-exclusive rule-based styling, 2 doesn't work: the feature counts in the legend show [0] for each AttType.

I can think of two workarounds, but neither is satisfactory:

One is to have 16 copies of the point layer, each filtered to a specific floor and named "Level 1 Points", etc., then setting up each map composer with a map ID as I already do ("Level 1"), only with a layer filter using a data-defined override (e.g., concat(@map_id,' Points')|'Other vector layer'|'Background raster')

However, having 16 point layers in the project makes it quite sluggish, and having to update any style changes in one layer across all 16 layers is something I'd really rather avoid.

The other is to manually generate feature counts in text boxes using aggregate() (e.g. aggregate('layer name','count',"AttType","AttType" = 'A' AND "Floor" = 1) and put them in the appropriate place next to each AttType on the legend.
This is even more time-consuming, but does avoid the problems of the latter workaround.

Is there any other way to solve this issue?
The Map Composer legend seems hard to customise.
My understanding is that I cannot use the Atlas function as the points all fall within the same XY plane (=building footprint).


texts: Identify identical geometry in column with common number in QGIS?

texts:
I'm trying to identify identical (or duplicate) geometry in QGIS by adding a common number for each identical set. For example, if there are 4 features, and 2 are identical:

1| Identical with 2

2| Identical with 1

3|

4|

I want the following:

1|1

2|1

3|2

4|3

With ArcGIS, there is a tool called "Find Identical" that creates a table for geometry duplicates - you choose a shapefile and choose "Shape" under the list to identify within, and it produces a table of identicals, as above. Here is a screenshot of the tool:


texts: how to use session variable in javascript to get parameter value from one javascript page to another?

texts:
I am trying to get value from one js file to another i am able to do that using cookies but i need to do using sessions can anyone suggest me on this ..(this value i need is for login control)
here is my code..

*********Login.js********

$(document).ready(function(){
    var assemblyCode;
    $('#btnSubmit').click(function () {
        var name = document.getElementById("txtUserName").value;

       // alert("cons:" + name);

        assmbylogin(name);

    });

         function assmbylogin(AssemblyName) {
           //  alert('hi');

             $.ajax({
                 type: "POST",
                 url: "Varuna.asmx/assmbylogin",
                 data: "{'AssemblyName':' " + AssemblyName + "'}",
                 contentType: "application/json; charset=utf-8",
                 dataType: "json",
                 async: false,
                 success: OnSuccess
             });

             function OnSuccess(data, status) {
               //  alert('hi');
                 var obj = JSON.parse(data.d);
                 var  assname, assemblynameArray = [];
                 for (var i = 0; i < obj.length; i++) {

                     assemblyCode = obj[i].AssemblyCode;
                     assname = obj[i].AssemblyName;
                     assemblynameArray.push(obj[i].AssemblyName);
                 }
                 alert(assemblyCode);
                 document.cookie = "" + assemblyCode + "; expires=Thu, 20 Aug 2016 12:00:00 UTC";

             }
  }
   });

*************health.js***********

$(document).ready(function () {
    var bloodsettings, pvthospsettings, diasettings, hospsettings;
    var bloodsettingsAssembly, hospsettingsAssembly, diapsettingsAssembly, pvthospsettingsAssembly;
    var statediacount, statebbcount, statephcount;

    var aseemblyCode = document.cookie;

    $("#health").hide();

    //health
    $("#hlt").click(function () {

        $("#bcHostel").hide();
        $("#education").hide();
        $('#groundwater').hide();
        $("#health").show();
        $("#barchart1").show();
        $("#barchart2").show();
        $("#barchart3").show();
        $("#barchart4").show();
        GetAssemblyDetails1();
        hospitalcategory();
        Hospitalcount();
        bloodbankscount();
        pvthospitalcount();
        diagnosticcenterscount();
       alert(aseemblyCode);

    });
});

in the above code am trying to get the assemblycode from login.js to health.js
(arcgis js api 3.17)


texts: Working with GeoServer and OpenLayers in Android?

texts:
I have an application using GeoServer and OpenLayers web platform, but must use the Android OpenLayers offline.

My doubt is the following, I can install GeoServer on Android? to work with OpenLayers using phonegap?

Do you have any suggestions on the other hand is the platform?

Dmitry Baryshnikov and west Efan provided links and explanations, most like to take one more question.

Today I have an application OpenLayers to be communicating with a server "GeoServer"

I would like to take my OpenLayers JavaScript in mobile, offline however, I would like somehow I store GeoServer data (WMS with JSON) that OpenLayers need on Android, or some other way to do this?


texts: Error in QGIS with Raster Tool in plugin

texts:
When I open QGIS I get the following error:

Couldn't load plugin 'rastertool' due to an error when calling its classFactory() method and the following command:

ModuleNotFoundError: No module named 'pyrsgis'
Traceback (most recent call last):
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 334, in _startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/lambr/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\rastertool\__init__.py", line 35, in classFactory
    from .RasterTool import RasterTool
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 744, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/Users/lambr/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\rastertool\RasterTool.py", line 54, in
    from pyrsgis.convert import rastertocsv
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 744, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'pyrsgis'

Python version: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS version: 3.10.1-A Coruña A Coruña, ef24c526da

Python Path:
C:/Users/lambr/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\leastcostpath
C:/OSGEO4~1/apps/qgis/./python
C:/Users/lambr/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/lambr/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
C:/OSGEO4~1/apps/qgis/./python/plugins
C:\OSGEO4~1\apps\Python37
C:\OSGEO4~1\apps\Python37\Scripts
C:\OSGeo4W64\bin\python37.zip
C:\OSGEO4~1\apps\Python37\DLLs
C:\OSGEO4~1\apps\Python37\lib
C:\OSGeo4W64\bin
C:\OSGEO4~1\apps\Python37\lib\site-packages
C:\OSGEO4~1\apps\Python37\lib\site-packages\win32
C:\OSGEO4~1\apps\Python37\lib\site-packages\win32\lib
C:\OSGEO4~1\apps\Python37\lib\site-packages\Pythonwin
C:/Users/lambr/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/lambr/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\aderyn
C:/Users/lambr/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\aderyn

I have seen that a lot of people have faced similar problems.
I downloaded OSGeo4W Network installer 64 bit. The last version of QGIS.

How can I solve this?

I ran the command but I got the following error:

Do you know something about it?

@SaijinNaib's answer worked, however when I open QGIS I got another error:

Couldn't load plugin 'rastertool' due to an error when calling its classFactory() method

ModuleNotFoundError: No module named 'visvis'

Do you know how to solve this as well?

Also why do I have so many errors?


texts: Warping a GeoTIFF File (Error1: Cannot open CSV)

texts:
Here is my code for the script :

import os
import subprocess as sub
# get your files
path = r'C:\Users\User\Downloads\jalna\pdf\test'
files = [os.path.join(path, f) for f in os.listdir(path) if f.endswith('.tif')]
# now run gdalwarp on each file
i = 0
for f in files:
    i+=1
    basename = os.path.basename(f)
    print(basename)
    newname = os.path.join('C:\Users\UserDownloads\jalna', str(i) + '_warped.tif')
    # set up call
    my_call = [r'C:\Program Files\QGIS 2.18\OSGeo4W.bat',
                'gdalwarp', f, newname, '-cutline', 'neat.csv', '-cwhere', 'id=1', '-crop_to_cutline']
    # call it
    print 'Using gdalwarp on {0} to convert it to {1}'.format(f, newname)
    p = sub.Popen(my_call, stdout=sub.PIPE, stdin=sub.PIPE)

    stdout, stdin = p.communicate()
    if p.returncode != 0:
        print stdout
        print stderr

I want to run this command from OSGeo4W Shell:

gdalwarp 1.tif 1_c.tif -cutline neat.csv -cwhere id=1 -crop_to_cutline

But getting the error:

  ERROR 1: Cannot open neat.csv.


texts: Does Google Earth Engine have InSAR Processing Capabilities for Sentinel-1 data?

texts:
I am attempting to find an easier way to work with InSAR data (i.e. getting land movement data). I have worked with SARProz and SNAP in the past, but was recently made aware that GEE had Sentinel-1 data. Through investigation, I found that GEE allowed crop monitoring and other types of monitoring that relied on backscatter intensities, but did not find any information about ground movement tracking.

Does GEE allow for ground movement (InSAR) work?


texts: Converting lat long strings to xy linestrings

texts:
In PostGIS, I have a dataset with road segments including a column with coordinates for the various segments. The coordinates are in a lat-long format and the data type of the column is character varying, plus it contains a varying number of coordinates in each row (see attached image).

What I wish to do is to be able to show the segments as geometric linestrings in order to later match them to a road network and/or a map.

I've already tried using ST_GeomFromText (see below) and ST_DumpPoint, but without any luck.

SELECT ST_GeomFromText('LINESTRING('||coordinates||')') as geom FROM vejnet.segmentnetdatex2format;

Which results in the following error code

  ERROR:  parse error - invalid geometry HINT:  "LINESTRING(55.67258,1"
  <-- parse error at position 21 within geometry SQL state: XX000

I'm still fairly new to PostGIS and working with GIS/spatial data in general.


texts: Adding custom Mapbox layer to Quickmapservices in QGIS

texts:
I'm trying to add a custom basemap, created in Mapbox, to Quickmapservices in QGIS.  However, the layer won't load and it keeps showing an error message:

Error: Layer xx can't be added to the map!

Filling out the Quickmapservices new WMS service form I've added:
-The URL for the WMTS link to Quickmapservices found at the bottom of 'share map' in Mapbox.
- The layer name (rather than the layer title), which I found by adding the basemap to QGIS through 'add WMS layer' first.

Does anyone know what might be going wrong? I've checked the QMS FAQ and haven't found an answer. Is it that a WMTS link won't work with WMS? Do I need to add any formatting to the link?


texts: Quantifying error on georeferenced aerial imagery

texts:
I georeferenced historical aerial imagery of a coastal area to estimate the size of oil sheens on the surface of the ocean. The control points I used were on land (street intersections, buildings, etc.), and I used thin plate spline interpolation. I don't know much about it besides what I learned from this page.

Is there any way I can quantify the error on the size of the oil sheens?


texts: arcgis server 10, vb.net add reference problems

texts:
i am using 'arcgis server 10' and 'visual web developer 2010 express'. now, i opened a 'web mapping application' then wanted to 'add a reference' to 'esri arcmap object library 10.0' it says:

 A reference to 'ESRI ArcMap Object Library 10.0' could not be added.
 Converting the type library to a .NET assembly failed. No process is
 associated with this object.

it said the same for 'esri 3danalyst object library 10.0', 'esri geoanalyst object library 10.0', 'esri geostatisticalanalyst object library 10.0', 'esri spatialanalyst object library 10.0'...

but for this one it says:

 A reference to 'ESRI Display Object Library 10.0' could not be added.
 Converting the type library to a .NET assembly failed. The type library
 importer encountered an error during type verification. Try importing
 without class members.

i tried this with target .net frameworks 3.5 and 4.0. plz can you explain how can i add these references successfully? thanks!


texts: spatialite.dll with Python and QGIS - the specified module could not be found

texts:
I'm having a bit of a nightmare getting to the bottom of this DLL loading problem. The stand alone script works absolutely fine but the exact same code does not work within QGIS.

import sys
from pysqlite2 import dbapi2 as sqlite

sys.path.append(os.path.dirname(__file__)
conn = sqlite.connect(db_name)
conn.enable_load_extension(True)
spatial_dll = os.path.join(os.path.dirname(__file__), "libspatialite-1.dll")
conn.execure('SELECT load_extension("%s")' % spatial_dll)

In the same directory as this script are the dependant DLL needed. It works like a charm stand alone, but not when put inside of QGIS. I am guessing I either have some kind of path issue or a compatibility issue with another version of the DLL somewhere else in QGIS.

Can anyone point me in the right direction? I don't have enough hair left as it is and I am losing more and more each hour...

Dan


texts: Can we count the points in polygon using pyqgis code?

texts:
Does anyone know how to calculate the point in polygon from the given points and the buffer output using python qgis.
I have every thing loading the data and retrieving the output. But now I want to count the points in polygon. I was checking this link but that didn't work.

Code:

    from PyQt4.QtCore import *
    from PyQt4.QtGui import *
    QgsApplication.initQgis()
    uri = QgsDataSourceURI()#4326

    #loading fisrt layer with multistring
    uri.setConnection("localhost","5432","db_name","postgres","password")
    uri.setDataSource("schema_name","table_name","column","")
    vlayer = QgsVectorLayer(uri.uri(),"LayerName1","postgres")
    QgsMapLayerRegistry.instance().addMapLayer(vlayer)

    #loading second layer with points around it

    uri.setConnection("localhost","5432","db_name","postgres","password")
    uri.setDataSource("schema_name","table_name","column","")
    vlayer = QgsVectorLayer(uri.uri(),"LayerName","postgres")
    QgsMapLayerRegistry.instance().addMapLayer(vlayer)

    #creating buffer
    mc = iface.mapCanvas()
    layer = mc.currentLayer()
    print layer.name()
    QgsGeometryAnalyzer().buffer(layer, "C:/Users/***/aa.shp", 0.0004, False, False, -1)

    #loading back buffer
    layer1 = iface.addVectorLayer("C:/Users/***/aa.shp", "layer_name", "ogr")
    layer1 = qgis.utils.iface.activeLayer()

#point in polygon
    provider = layer.dataProvider()
    provider.select()
    feat = QgsFeature()
    vertex_count = 0
    while provider.nextFeature(feat):
          geom = feat.geometry()
          if geom.type() == QGis.Polygon:
            layer_vertices = 0
            if geom.isMultipart():
              polygons = geom.asMultiPolygon()
            else:
              polygons = [ geom.asPolygon() ]
            for polygon in polygons:
              for ring in polygon:
                layer_vertices += len(ring)
            print "Feature %d: Vertices %d" % ( feat.id(), layer_vertices)
            vertex_count += layer_vertices
            print "Total vertices: %d" % (vertex_count)

ERROR:

provider.select()
AttributeError: 'QgsVectorDataProvider' object has no attribute 'select'

Also please can you guys share the link or the code for counting points in polygon. Few of them posted the same error can anyone share the code with explanation for point in polygon.


texts: Importing a shapefile into Postgres: Can I use Multiple CPU Cores?

texts: I'm importing some shapefiles into Postgres/PostGIS.  I'm using Windows 8 and the shapefile importer/exporter GUI.  It maxes out one CPU core. I'm wondering if there is a way I can use multiple cores to increase the speed?


texts: GetFeatureInfo for gxp_olsource in OGS GeoExplorer

texts:
We have an OpenGeo Suite GeoExplorer 3.0.1 and want to make GetFeatureInfo-requests, via the identify button, for WMS-layers of type gxp_olsource but it doesn't work, I get no response at all and I don't get any errors in the browser debug console. It works great to make a GetFeatureInfo-request for WMS-layers with sources of type gxp_wmscsource
or gxp_wmssource.

Is it supposed to be working in the original version of GeoExplorer or do I have to make adjustments to the code?

Layer configuration:

...
sources: {
...
           ol: { ptype: "gx_olsource" }
},
map: {
...
        { group: "groupfisketsgeografier",
        visibility: false,
        queryable: true,
        source: "ol",
        type: "OpenLayers.Layer.WMS",
        args: ["Fredningsomraden",
            "http://geodata.havochvatten.se/geoservices/hav-fisketsgeografier/ows",
            { layers: "hav-fisketsgeografier:havet-fredningsomraden",
            transparent: true,
            isBaseLayer: false,
            tiled: true} ]
            },

        }
...


texts: Creating Index Map for ArcGIS Pro Map Series?

texts:
I have used the steps below to try and create a nice looking Index Map for a test Map Series in ArcGIS Pro 1.4.1 using an index which is a regular grid.  The idea is to have one map page in the large map frame and then to have an index in a small map frame that keeps the main map page centred and shows the eight surrounding map pages (with them labelled).

Start ArcGIS Pro and create a new Map project as C:\Temp\TestProj
Use the Python pane to run the four lines of code below to create a fishnet called TestFishnet in C:\Temp\TestProj\TestFishnet, give it a Geographic Coordinate System (GDA94) and add/calculate a field called PageName:

arcpy.management.CreateFishnet(r"C:\Temp\TestProj\TestProj.gdb\TestFishnet", "110 -35", "110 -20", 2, 1, None, None, "150 -20", "NO_LABELS", "DEFAULT", "POLYGON")
arcpy.management.DefineProjection(r"C:\Temp\TestProj\TestProj.gdb\TestFishnet", "GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
arcpy.management.AddField(r"C:\Temp\TestProj\TestProj.gdb\TestFishnet", "PageName", "TEXT", None, None, 10, None, "NULLABLE", "NON_REQUIRED", None)
arcpy.management.CalculateField(r"C:\Temp\TestProj\TestProj.gdb\TestFishnet", "PageName", '"Page"+str(!OID!)', "PYTHON_9.3", None)

Zoom To Layer on TestFishnet and make its symbol a Black Outline (1pt)
INSERT | New Layout > ISO - Landscape A3
INSERT | Map Frame > Map Default
LAYOUT | Map Series > Spatial and accept defaults (Map Map Frame, TestFishnet, PageName, PageName) then OK
Change the Map Frame on the Layout so that it hugs the boundary of the map sheet and move to leave room for an Index Map (e.g. for Page100 see image below)

To prepare an index map frame start with INSERT | New Basemap
Drag and drop C:\Temp\TestProj\TestProj.gdb\TestFishnet onto the Basemap and Label that layer with the PageName
On the Layout, INSERT | Map Frame > Basemap Default and move it to get something like the image below

To configure the index map frame, in the Contents pane use List By Map Frame to right-click and open Properties of Basemap Map Frame
On the Display Options tab set Constraint to Linked map frame extent and Map Frame to Map Map Frame, leave Layer as  and set Margin to 200 Percent.
By manually resizing the index map frame I can almost get what I want (below) but the labeling is very messy.

Is it possible to improve on this so that in the index map frame it shows 9 labelled map sheets that update nicely as I step through the pages of the map series?


texts: AttributeError: 'NoneType' object has no attribute 'length' when running function for compactness

texts:
I want to calculate compactness indices for geometries I have in a geodataframe.
I'm using the following function  as was mentioned here (Measuring compactness in Python)  by Kadir in order to calculate different shapes compactness. This is the function:
def pp_compactness(geom): # Polsby-Popper
    p = geom.length
    a = geom.area
    return (4*pi*a)/(p*p)

def s_compactness(geom): # Schwartzberg
    p = geom.length
    a = geom.area
    return 1/(p/(2*pi*sqrt(a/pi)))

df['polsby_popper']=df['geometry'].apply(pp_compactness)
df['Schwartzberg'] = df['geometry'].apply(s_compactness)

This is my dataframe sample:
>>>ID  geometry                  centroid              class        perimeter
0  43 POLYGON((9.956 9.12...)   POINT (9.95 9.10)      Forest       521.347
1  46 POLYGON((8.941 8.54...)   POINT (8.5 9.44)       Desert       345.123
...

I have tried to apply the function on the geometry column as following:

f['polsby_popper']=df['geometry'].apply(pp_compactness)
df['Schwartzberg'] = df['geometry'].apply(s_compactness)

but I'm getting the following error message:

AttributeError: 'NoneType' object has no attribute 'length'

I have tried to set the geoemtry before and run again:
df=df.set_geometry('geometry')

but still got the same error message.
I have no idea why it is missing the length attribute,
My end goal is to use this function to calculate compactness in geometries I have as GeoPandas .
EDIT: I have changed the function instead of geom.length and geom.area to geometry.length and geometry.area (because that's the name  of the geometry column) but that raised the same error.
def pp_compactness(geometry): # Polsby-Popper
    p = geometry.length
    a = geometry.area
    return (4*pi*a)/(p*p)

def s_compactness(geometry): # Schwartzberg
    p = geometry.length
    a = geometry.area
    return 1/(p/(2*pi*sqrt(a/pi)))


texts: QGIS, WKT (Well Known Text) LINES and POLYGON files loaded ok, but POINT not

texts:
I have several files in WKT format, and something strange happens when I import them in QGIS (v.2.18.1) (layer->add layer->add delimited Text Layer), files containing lines and polygons are correctly located on the map, but those containing points are not, and the same CRS is used  (EPSG: 4326, WGS84) for all files.

The examples that I put correspond to data an area of the north of Spain, and as you can verify, the lines and points are loaded in that zone, but the points go to the Indian Ocean, in the coast of Tanzania.

Here you can download the wkt files used (are only 65KB zip file) : https://drive.google.com/open?id=0B2tCGwcrrXRtN2tHSldfbDJ6MlE

It is very strange, because the coordinates are good, if I take a point of the file that is misplaced, and I copy the coordinates within google maps, the point is correctly located in the north of Spain.

Here I put some lines (only the first ones) with the content of these files (utf-8).

File-> Point.wkt

type|label|geometry
0x6510|Océano Atlántico|POINT(43.71632 -8.14070)
0x6510|Océano Atlántico|POINT(43.73159 -7.99925)
0x6510|Océano Atlántico|POINT(43.80043 -7.98792)
0x5000||POINT(43.66791 -8.05667)
0x5000||POINT(43.67915 -8.03564)
0x5000||POINT(43.68267 -7.99418)
0x6400|Parque Eólico Couce Do Penido|POINT(43.67048 -7.95805)
0x6400|Parque Eólico Couce Do Penido|POINT(43.67100 -7.93513)
0x5000||POINT(43.70138 -8.01916)

File-> Lines.wkt

type|label|geometry
0x15||LINESTRING(43.69263,-8.05847 43.68713,-8.06396 43.68164,-8.05847 43.68164,-8.06946 43.67615,-8.06396)
0x15||LINESTRING(43.67615,-8.06396 43.67615,-8.06946 43.67615,-8.06396 43.67615,-8.06946 43.67065,-8.06946 43.66516,-8.06946)
0x14||LINESTRING(43.66516,-7.85522 43.67065,-7.84973)
0x15||LINESTRING(43.70911,-8.05298 43.70361,-8.05847 43.69812,-8.05847 43.69812,-8.06396 43.69263,-8.05847)
0x15||LINESTRING(43.71460,-8.03650 43.70911,-8.04199 43.71460,-8.04199 43.70911,-8.04199 43.71460,-8.04749 43.70911,-8.05298)
0x15||LINESTRING(43.71460,-8.01453 43.70911,-8.02002 43.71460,-8.02551 43.71460,-8.03650)
0x15||LINESTRING(43.70911,-8.00903 43.71460,-8.01453)
0x15||LINESTRING(43.71460,-7.97058 43.70911,-7.98706 43.70361,-7.99805 43.70361,-8.00354 43.70911,-8.00903)
0x15||LINESTRING(43.70911,-7.85522 43.70911,-7.84973)

and File-> Polygons.wkt

type|label|geometry
0x32||POLYGON((43.69263,-8.06946 43.69263,-8.06396 43.68164,-8.06396 43.68164,-8.06946 43.67065,-8.06946 43.66516,-8.06946 43.66516,-8.18481 43.69263,-8.18481))
0x32||POLYGON((43.68164,-8.06946 43.67615,-8.06396 43.67615,-8.06946 43.67065,-8.06396 43.67615,-8.06946))
0x32||POLYGON((43.69263,-8.06396 43.69263,-8.05847 43.68713,-8.06396 43.68164,-8.06396))
0x4b|0001 Cariño|POLYGON((43.74756,-8.18481 43.74756,-8.01453 43.66516,-8.01453 43.66516,-8.18481))
0x32||POLYGON((43.71460,-8.04199 43.70911,-8.05298 43.71460,-8.05298))
0x32||POLYGON((43.70911,-8.00354 43.71460,-8.01453 43.70911,-8.02002 43.71460,-8.02551 43.71460,-8.00354))
0x32||POLYGON((43.71460,-8.00354 43.71460,-7.97607 43.70361,-7.99805 43.70911,-8.00354))
0x32||POLYGON((43.70911,-7.85522 43.70911,-7.84973 43.70361,-7.84973))
0x4b|0001 Cariño|POLYGON((43.74756,-8.01453 43.74756,-7.84973 43.66516,-7.84973 43.66516,-8.01453))


texts: Using ArcGIS field calculator with if else to populate field?

texts:
ArcMap 10.4.1 on Windows 7.

I am trying to fill in some blank rows in a column using field calculator. The field I am trying to fill is: label. If it is empty I want it to be populated with a value from the column: FULLNAME. If it is not blank I want it to retain its current value.

I have looked at other if else examples and this seems to be very similar, but my code doesn't work. I am also confused on when the exclamation points, example (!field!) are necessary if field calculator expressions.

Do those tell python to pull that field from the table?

def updateLabel (label, FULLNAME):
  if !label! is None:
     return !label!.replace( !FULLNAME! )

Do I need a code block in this case?


texts: Openlayers 3 equivalent of render intents of OpenLayers.Stylemap?

texts:
I want to convert my openlayers 2 features to openlayers 3. I have following part in openlayers 2. I need an equivalent of StyleMap in openlayers 3 so that I can distinguish the default action and select action.

var subMap = new OpenLayers.StyleMap({
    "default": new OpenLayers.Style({
        fillColor:"orange",
        strokeColor:"orange"
    }),
    "select": new OpenLayers.Style({
        fillColor:"#00FFFF",
        strokeColor:"#00FFFF"
    })
});

I can write the following for "default" case. But where to mention for the "select" case?

var subMap= new ol.style.Style({
    fill: new ol.style.Fill({
        color: "orange"
    }),
    stroke: new ol.style.Stroke({
        color:"orange"
    })
});

Can anyone please help in this regard?


texts: QGIS delimited text import CRS not being defined

texts:
I've recently updated to both the long term release and 2.10 versions of QGIS and the delimited text plugin has stopped working.

Everytime I try to add a CSV file and say ok on the plugin it no longer asks for me to define the CRS and then defaults to WGS 84 (image below). Given I work in UTM it's a massive hassle to convert to long/lat everytime just so it wil display correctly so i can transform it to UTM.

Any suggestions?


texts: Changing field type of Excel attribute table for a Join

texts:
I have a string field type of shapefile already in ArcMap. I need to import an excel file in and join with the shapefile. The only problem is the field type of the shapefile is a String and the field type for the excel file is a Double.

The excel field is all numbers, but the shapefile field is mostly numbers and a few letters mixed in, I think thats why it came out as a string.

I know you can't add a field or use the field calculator into a table exported from excel, I'm currently using it as a .csv.


texts:
I have a geopackage for India taken from GADM. I can plot all the districts of India by using layer 1, all states using layer 3, and the whole country using layer 4. Here is the code:

> library(rgdal)
Loading required package: sp
rgdal: version: 1.3-3, (SVN revision 759)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: C:/Users/H.P/Documents/R/win-library/3.4/rgdal/gdal
 GDAL binary built with GEOS: TRUE
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: C:/Users/H.P/Documents/R/win-library/3.4/rgdal/proj
 Linking to sp version: 1.3-1
Warning message:
package ‘rgdal’ was built under R version 3.4.4
> layers<-ogrListLayers("gdam36_IND.gpkg")
Error in ogrListLayers("gdam36_IND.gpkg") : Cannot open data source
> layers<-ogrListLayers("gadm36_IND.gpkg")
> x<-readOGR("gadm36_IND.gpkg",layers[3])
OGR data source with driver: GPKG
Source: "C:\Users\H.P\Documents\gadm36_IND.gpkg", layer: "gadm36_IND_1"
with 36 features
It has 10 fields
> plot(x)

Which outputs :

I can get the summary of the variable x by:

> summary(x)
Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x 68.186249 97.41516
y  6.754256 35.50133
Is projected: FALSE
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
 GID_0      NAME_0        GID_1                    NAME_1
 IND:36   India:36   IND.1_1 : 1   Andaman and Nicobar: 1
                     IND.10_1: 1   Andhra Pradesh     : 1
                     IND.11_1: 1   Arunachal Pradesh  : 1
                     IND.12_1: 1   Assam              : 1
                     IND.13_1: 1   Bihar              : 1
                     IND.14_1: 1   Chandigarh         : 1
                     (Other) :30   (Other)            :30
                                                                                VARNAME_1
 Ãles Laquedives|Laccadive|Minicoy and Amindivi Islands|Laccadives|Lackadiverna|Lakkadiven|Lakkadi  : 1
 Agence de la FrontiÃ¨re du Nord-Est(French-obsolete)|North East Frontier Agency                     : 1
 Andaman & Nicobar Islands|Andaman et Nicobar|Iihas de Andama e Nicobar|Inseln Andamanen und Nikobare: 1
 Bangla|Bengala Occidentale|Bengala Ocidental|Bengale occidental                                     : 1
 DAdra et Nagar Haveli|Dadra e Nagar Haveli                                                          : 1
 (Other)                                                                                : 9
 NA's                                                                                :22
 NL_NAME_1            TYPE_1             ENGTYPE_1    CC_1        HASC_1
 NA's:36   State         :28   State          :28   NA's:36   IN.AN  : 1
           Union Territor: 8   Union Territory: 8             IN.AP  : 1
                                                              IN.AR  : 1
                                                              IN.AS  : 1
                                                              IN.BR  : 1
                                                              IN.CH  : 1
                                                              (Other):30

Is there a way to plot only for a particular state? That is, for example, x$NAME_1=="Chandigarh"?


texts: JSONP working with GeoServer and OpenLayers

texts:
I'm in this scenario (fighting with the same origin police issue):

I have a map initialized with OpenLayers and a WFS request trough the Protocol.Script object that is giving me this error:

Uncaught TypeError: Object function (a,b){this.metadata={};b=OpenLayers.Util.extend({},b);null!=this.alwaysInRange&&(b.alwaysInRange=this.alwaysInRange);this.addOptions(b);this.name=a;if(null==this.id&&(this.id=OpenLayers.Util.createUniqueID(this.CLASS_NAME+"_"),this.div=OpenLayers.Util.createDiv(this.id),this.div.style.width="100%",this.div.style.height="100%",this.div.dir="ltr",this.events=new OpenLayers.Events(this,
this.div),this.eventListeners instanceof Object))this.events.on(this.eventListeners)} has no method 'addOptions'

The error seems to be related with OpenLayers lib (maybe the version, actually is the 2.13), but i don't know.

What i tried is:

consulados = OpenLayers.Layer.Vector("Consulados", {
        strategies: [new OpenLayers.Strategy.BBOX()],
        protocol: new OpenLayers.Protocol.Script({
            url: "http://"+ip+':'+port+"/geoserver/wfs",
            callback: 'handleResponse',
            params: {
                service: "WFS",
                version: "1.1.0",
                srsName: "EPSG:4326",
                request: "GetFeature",
                typeName: "topp:tazmania_roads",
                outputFormat: "text/javascript"
            },
        })
});


texts:
I downloaded OSM data for Austria from Geofabrik and converted them with ogr2ogr into a SpatiaLite database so that I could use some Google Maps style sheets I found.

Until here, everything worked very well.

Then I wanted to add a shapefile with some points I tracked in nature as well as a orthophotografic I bought at an official government store.

That works also, but the OSM data on the one hand and everything else on the other do not match. My own data is about 50 meters farther east than it should be to match the correct point on the OSM data.

I looked for on-the-fly transformation, but as I read this feature isn't optional any more since version 3.

I use QGIS 3.0.1 and as CRS WGS 84 (EPSG 4326)


texts: How can I disable file path balloon in Layer panel? When I pass the cursor in any object in Layer panel I see an annoying long text showing the file path. I would like to disable it.


texts:
I wan to get 256px tiles of a geoserver wms service. I specified a bounded box like this:

-87.746754646322, 41.894876003264 -87.710705757162, 41.869813442229

But tiles are pyramid format. Wms service is 9 level and request parameter does not contain zoom level parameter.

http://localhost:8080/geoserver/GIS/wms?LAYERS=GIS:ways&STYLES=&FORMAT=image%2Fpng&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG%3A36100&BBOX=-87.746754646322,41.894876003264-87.710705757162,41.869813442229&WIDTH=259&HEIGHT=256

This request returns all ways that bbox contains.

Bu I want to get all tiles with programatically like pyramid in bbox like this:

tile1.png
tile1.png tile2.png
tile1.png tile2.png tile3.png
tile1.png tile2.png tile3.png tile4.png


texts: On a pedestrian wayfinding scenario, I need to drape some extruded buildings close to the road as the user navigates along the road on a 2D Map Background. I hope I can use OpenLayers and ReadyMap for this. I have vector data in shape file format. I am a newbie to web-mapping. Any ideas on this to start with?


texts:
I have a shape with 110 polygons (hydrographic basins in certain computational sections).  I have to sum up their surfaces from upstream to downstream. below is an example of polygons representation. for example, sum for polygon 2 is  surface 1 + surface 2, for polygon 4 is surface 1 + surface 2 + surface 3 + surface 4.... and so on.

Is there any Python script or a ModelBuilder to automatically determine the surface for all polygons in the computational sections?


texts: GEE filterBounds not working with S5P dataset

texts:
It seems to me ee.ImageCollection.filterBounds is not filtering properly a dataset from Sentinel-5P.
I get images which do not include the point, which is the same for filterBounds and Map.setCenter.

var dataset = 'COPERNICUS/S5P/OFFL/L3_NO2';
var band = 'tropospheric_NO2_column_number_density';
var begin_date = '2019-05-03'; var end_date = '2019-05-04';
var collection = ee.ImageCollection(dataset)
  .filterDate(begin_date, end_date);

// Number of S5P NO2 images in one day
print(collection.size()) // 14

var lat=18.232527; var long=-66.257565; // Puerto Rico
var point1 = ee.Geometry.Point([long, lat])
var collection_clean = collection.filterBounds(point1)

// Number of S5P NO2 images in one day
print(collection_clean.size()) // 14????? All images cross Puerto Rico????

// Take the first image and try to display it
var image = collection_clean.first().select(band);
var zoom=8;
Map.setCenter(long, lat, zoom);
var min_value = 0.00001; var max_value = 0.00003;
var visParams = {min: min_value,
  max: max_value,
  palette: ['green', 'blue', 'yellow', 'red']}

// Image can only be seen after zooming out 5 times
Map.addLayer(image, visParams, "NO2")

If you zoom out 5-6 times you can finally see the NO2 layer.


texts: Is there any way to find the shortest distance from a point to a line?

texts:
I´m working with QGIS 1.8.0. I have a map divided into many features (thousands of small rectangles) and I have a point at the center of each of those features. I also have a map of the transmission lines of the country. Is there any method to measure (not one by one) the shortest distance from the points to the transmission lines? Does the maps have to be projected into any specific CRS?

Regards


texts: Select by location not selecting all features in ArcPy

texts:
I am having an issue with the "select layer by location" tool in my ArcPy script, that may be due to a hidden memory error. Using one feature class with 281 polygons, I would like to select all of the possible bird ranges that overlap each polygon. There are ~12000 bird range polygons. My goal is a table with all of the 281 polygons and all the bird ranges associated with each one. Simply running the select by location tool runs into a memory error so I am attempting to iterate through each of the 281 shapes and export them to text file to later merge together.

The first 4 polygons accurately capture all of the intersecting polygons. Oddly, after that the script does not crash, but leaves out several hundred bird ranges for each of the 281 polygons. I know this by using the select by location tool on one feature in ArcMap and comparing the results.

Any ideas?

Maybe there is a different approach I haven't thought of.

def makeFeatureLayer():
    env.workspace = pathlist[0] # Set initial workplace to file geodatabase
    BirdData = "Joined_Study_Data"
    tempBirdData = "BirdsDataTemp"
    #selection is performed on a temporary layer
    arcpy.MakeFeatureLayer_management(BirdData, tempBirdData)
    print "Made temp layer:",tempBirdData
    return tempBirdData

def selectionToText(temp_layer):
    env.workspace = pathlist[0] # Set initial workplace to file geodatabase
    FPUData = "FPU_proj"
    cursor = arcpy.SearchCursor(FPUData)
    for count,row in enumerate(cursor): # loop on each FPU feature
        temp_layer = makeFeatureLayer() # make temp layer to avoid mem error
        arcpy.SelectLayerByLocation_management (temp_layer, "INTERSECT",
                                                #select based on geometry
                                                row.SHAPE)
        outname = pathlist[1]+"JoinSelection"+str(count)+".txt"
        #export selection to text file, specifying fields of interest
        arcpy.ExportXYv_stats(temp_layer,["SPCRECID","PRESENCE","ORIGIN",
                                            "SEASONAL", "scientific_name",
                                            "common_name", "rl_category"],
                                "COMMA", outname, "ADD_FIELD_NAMES")
        print "Finished intersecting birds with FPU num:", str(count)
        count+=1

def main():
    temp_layer = makeFeatureLayer()
    selectionToText(temp_layer)


texts: How to render GeoJSON tiles from Tilestream for polymaps?

texts:
Have been using polymaps successfully till now.

Need to generate GeoJSON tiles for Canada Provinces.

Found Tilestream and got it working with Tilemill.

Is there a way I can make Tilestream hand out GeoJSON tiles for polymaps without Wax?

http://my_tilestream_server:8888/provinces/{Z}/{X}/{Y}.json

EDIT

Any help here would be greatly appreciated.

Still trying to make either Tilestream / Tilemill or Tilestache work correctly.

Does anyone know where I can find a good shapefile for Canada and its Provinces (Canadian Fips code would be a bonus)?  Only ones I have found are not accurate at zoom levels greater than 8.

1 - Tilestream / Tilemill are great, but need a way to dynamically add additional data to the Hover.  Can't seem to find a way to do this without including it in the dbf.

2 - Installed Tilestache, but having problems converting data from shapefile to sql.  I think I am messing missing something in the conversion to sql.  Or is there a data source I can import that already has this info?


texts: Converting 0 values to 1 in QGIS 3 using r.mapcalc GUI

texts:
I'm using GRASS 7.4.1 in QGIS 3.2 and would like to change 0 values of a slope raster to 1. I've been following this video tutorial and my r.mapcalc function looks exactly the same. However the output adds 1 to all values in the slope raster. The model looks like this in the GUI:

And the formula in the output reads "slope_1 = ((slope+((slope==0)+1)))"

How can I convert only 0 values to 1 instead of adding 1 to all values? In ArcGIS raster calculator I would use something like Con(raster==0,1,raster)

Stripping the GRASS r.mapcalc function back further the output from this created a raster with original 0 values converted to 2 and all other values converted to 1:

Update

I am aware there are other methods I could use to achive the same results but part of the purpose of doing this is to learn how to use the r.mapcalc GUI and so I am only looking for answers that use r.mapcalc GUI. Specifically I'd like to know why the formula / model used isn't behaving as expected and if another formula / model would be better suited to the problem.


texts: What are the best ways to connect two coordinates on mapinfo?

texts: I have with the help of this forum drawn a coordinate system with longitude and latitude at different location. I am trying to connect these sites what is the best way to connect. Is it just drawing lines or something else?


texts: Button on click does not display pie charts

texts:
I am trying to display a pie chart on the click of a button using JavaScript/HTML. The pie chart displays fine without button click. But it doesn't work with the button. I have tried the following:

<button class="btn btn-danger navbar-btn" onclick="drawChart(new google.visualization.PieChart(document.getElementById('chart_div')))">Button</button>

google.charts.load('current', {'packages':['corechart']});

      // Set a callback to run when the Google Visualization API is loaded.
google.charts.setOnLoadCallback(drawChart);
google.charts.setOnLoadCallback(drawAssistChart);

      // Callback that creates and populates a data table,
      // instantiates the pie chart, passes in the data and
      // draws it.
    function drawChart() {
    // Create the data table.
    var data = new google.visualization.DataTable();
    data.addColumn('string', 'Golazo');
    data.addColumn('number', 'Gol');
    data.addRows([
        ['Messi', 12],
        ['Eto,o', 8],
        ['Gomez', 8],
        ['Ronaldo', 6],
        ['Soldado', 6],
        ['Raul', 5]
        ]);
        // Set chart options
    var options = {'title':'Top Goal Scorers 2010-11',
    'width':300,
     'height':200};

// Instantiate and draw our chart, passing in some options.
//var chart = new google.visualization.PieChart(document.getElementById('chart_div'));
//chart.draw(data, options);
}

Could someone kindly suggest a solution?


texts: 2 Similar Queries Run at Vastly Different Speed in a Very Surprising Way

texts:
I have to be very careful here. Both queries are so similar that the only different is the area searched.

What's surprising is the query that run on BIGGER AREA, and naturally yield MORE ROWS, is the one that run way FASTER.

In fact the only reason why I want to run the second query is to narrow down the search because the original query returns too many results.

SELECT BusinessID
FROM tableauxiliary
WHERE MBRWithin( Latlong, GeomFromText( 'MULTIPOINT(-6.2209869856406 106.73701301436,-6.1490130143594 106.80898698564)' ) )
AND MATCH (
FullTextSearch
)
AGAINST (
'res*'
IN BOOLEAN
MODE
)

returns

  Showing rows 0 - 29 ( 5,007 total, Query took 0.0226 sec)

Now my code will streamline things up.

I asked the same query except for SMALLER regions

SELECT BusinessID
FROM
  tableauxiliary
WHERE
  MBRContains(
        GeomFromText (
            'MULTIPOINT(-6.1928749092968 106.7651250907,-6.1771250907032 106.7808749093)'
            ),
            Latlong)
  AND MATCH (FullTextSearch) AGAINST ('res*' IN BOOLEAN MODE)
AND Prominent >15

  Showing rows 0 - 29 ( 30 total, Query took 15.3889 sec)

What gives?


texts: Query M value where point intersects route in ArcGIS?

texts:
Is there a way I can query the M value along my route where a point touches the line. I'd like to use the field calculator in my Points fc to query the M value where the point intersects the line and populate the field with an M value. My points are M aware.

I've done something slightly similar by querying the beginning and end measure values from a line and applying it to a field in the line's table using !Shape.extent.MMin!.

Rather than manually inputting the measure value where the point and line intersect using the Route Identify tool, I'd like to automate this process if possible.


texts: How to add link in ArcGIS Online web map popup to google or bing driving directions?

texts:
I have an ArcGIS webmap with school boundaries. I would like to provide a link in the pop up that takes a public user to a free routing website like Bing or Google with the directions to their neighborhood school already filled in.

http://agtrail.rivcoca.org/ The County of Riverside has what I am talking about on their webmaps.

Is this something I will have to create a url for for each school and add it to each school boundary attribute? Or is there a faster way?

I know ESRI can do routing directions but they use up credits and I want to do this without costing my district any more money.

My scripting skills are very basic but if coding is the only way to do this I am willing to try my best.


texts: Merging Features with Equal Values in one Attribute Table Field?

texts:
I'm trying to merge all the features in a shapefile that have the same value under a certain field in the attribute table. Right now, I open the attribute table and use "Select by Attributes" to select all the features of a certain value, merge those, then move on to the next value.

Is there a way to do this automatically?


texts: Find points that lie within a set of coordinates

texts:
I have a PostGIS db, and I want to find the points that lie within some region (a bounding box). I have two sets of coordinates, from which I can get all four rectangle points that form the box (if needed). My data column in question is named 'point' and it is also of type point.

Is there anyway to specify four sets of coordinates (lat/long) and get all the points that lie within the box.
Or specify two points and let the DB work out the rectangle's corners, and return the points within

Just in case I am not being clear as to what I want to achieve. The equivalent with 'vanilla' sql if I had a lat and long field instead of a point would be:

SELECT * FROM myTable where lat> xMin AND lat < xMax AND long > yMin and long < yMax

UPDATED EDIT:

I am trying underdark's solution. At first I didn't have the ST_MakePoint constuctor (now I do) and I still get a very similar error (just on a different character).

SELECT * FROM myTable WHERE ST_Within(ST_MakePoint(point),GeometryFromText('POLYGON((75 20,80 30,90 22,85 10,75 20))',4326))

and I am getting this error:

ERROR:  function st_makepoint(point) does not exist
LINE 1: SELECT * FROM triples WHERE ST_Within(ST_MakePoint(point),Ge...
                                          ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

********** Error **********

ERROR: function st_makepoint(point) does not exist
SQL state: 42883
Hint: No function matches the given name and argument types. You might need to add explicit type casts.
Character: 39

EDIT:

In the short term I can solve it with:

SELECT * FROM triples WHERE box '((point1),(point2))' @> point

But I will have to work out why none of the PostGIS functions are not working for me.


texts: Where is date information in OpenStreetMap data files?

texts:
I would like to access temporal information from OSM datafiles (not the changeset file)
Files *.osm.pbf downloaded from the GeoFabrik website or another OSM mirror do not show the following information:

Creation date
Last update date

I treat the data using gdal ogr*
How can I retrieve this information ?


texts: Preventing Calculate Field tool slowing when iterating in ModelBuilder?

texts:
I am using ModelBuilder to get some results (using various tools) and assign the results to a .dbf file. Each time the For loop runs the program calculates some values and assigns it to designated rows of the export.dbf file. For that I am using "Calculate Field" Tool. This tool assigns the calculated value i.e OFIPS, DFIPS, SNO, NearFID and NearDist to designated rows of Export.dbf. The tool works really quick (1-2 secs) at first but keeps on slowing down (8-10 secs) after say, hundreds of iterations. The speed of the "Calculate field" tool is critical to the speed of the overall ModelBuilder.

Is there any way that the ModelBuilder not slow down the time for each iteration as the number of iterations increase.

The above screenshot shows the slow part in blue circle.


texts: Split lines and sequentially number each fragment in QGIS

texts:
I have a line layer representing watermains that I need to split long lines into sections.  Each line has a unique ID number, but upon splitting a line this attribute is no longer unique as the fragments each have the same attribute.  I would like to be able to number each of the unique fragments sequentially so that we can identify them as originally belonging to a larger parent, but still be unique.  Upon splitting, each line breaks into different number of smaller sections.
For example

Line 1 becomes lines 1.1, 1.2, 1.3
Line 2 becomes lines 2.1, 2.2, 2.3.....2.47
etc

I can see how to number the whole of the data set, but not number individual subsets.  Having the number as a new attribute would be acceptable.


texts: Setting the same min/max for multiple raster layers QGIS

texts:
I have 268 hgt-files from Viewfinder Panorama that I'm using to create an elevation model, but each of the files has its own min and max-value in Layer Properties - Symbology - Band Rendering.

Currently it looks like this:

I found I can set the min and max of each layer manually, to get an uniform scale, but this is a hassle for all 268 layers. So I was wondering if I'm missing an option somewhere, or if im doing something completely crazy?


texts: Animating alterations to DEM based on function using ArcGIS for Desktop?

texts:
I've generated a function that alters height values of a DEM based on certain criteria. What I'm looking at doing now is animating this in ArcScene. Is there a way, in ArcGIS, for each frame to be the result of a raster calculation?

eg Frame_Timecode-->(function)-->Adjusted_DEM-->render

I'm a strong Python programmer, but have not yet done any scripting with ArcGIS, so I'm quite happy if the answer involves some coding.


texts: How do I generate a sampling path by numbering points?

texts:
I need to number points from left to right and then the next row of points the numbering must continue from right to left, then again from left to right the row thereafter - and so on, and so on.

How do I do this in Qgis?  The points are not necessarily in straight lines.


texts:
I am trying to pull large time series of ECMWF data over single pixels. I'm finding this to be a very slow operation and my guess is that it has to do with the mismatch between how the data is stored natively on platform (WGS) vs. how tiles are defined (UTM), making the raster operation very slow.
I am defining a single-pixel tile as so:
tile = dl.scenes.DLTile.from_latlon(lat, lon, 30000, 1, 0)
Even when I try to pull a single day's worth of data (10 bands, 24 scenes in a day) it takes > 1 minute for the raster call to succeed. This won't scale as I need to pull years worth of data over thousands of these single pixel areas.
In an attempt to do everything in WGS, I tried something like:
aoi = dl.scenes.AOI(geometry=geom, crs='EPSG:4326', resolution=0.25)
where my geom is a GeoJSON Point geometry. However, this didn't work as no bounds are able to be computed for a single point.
Is there a better way to be pulling ECMWF data?


texts: PyQGIS add feature in layer with Value Relations in columns

texts:
Have a point layer (from geopackage) with value ralations widget in several columns. Trying to add new feature to this layer with the following code:
import uuid
layers = QgsProject.instance().mapLayersByName('L10_Public--Point')
layer = QgsVectorLayer(layers[0].dataProvider().dataSourceUri(), '', 'ogr')
print(layer)

caps = layer.dataProvider().capabilities()

if caps & QgsVectorDataProvider.AddFeatures:
    feat = QgsFeature(layer.fields())
    feat.setAttributes([2, 'Public', NULL, str(uuid.uuid4()), str(uuid.uuid4()), '602010904', NULL, NULL, '97629445', NULL, 5.0, 2.0, 3.0, 4.0, NULL, NULL, 4.0, NULL, NULL, NULL, NULL, 4.0, 4.0])
    feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(1605057.281,645118.237)))
    res, outFeats = layer.dataProvider().addFeatures([feat])

And getting:
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.10\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 2, in <module>
IndexError: list index out of range

Should be mentionel that the same code works well for another layer, which has no widgets in columns. Anybody knows what IndexError means here? Or how can I programmatically add new features in such a layer?


texts: Individual tree segmentation workflow with lidR package

texts:
I am seeking advice on how to build a workflow for individual tree segmentation (ITS) processing with the lidR package using catalog objects.

My current ITS workflow is parallel and leverages the lidR package, but it does not use LAScatalog objects. I am not quite sure how to make the LAScatalog object with a directory of DTMs and CSMs to use in a workflow like the one described here

All of the documentation for ITS in the wiki appears to assume that the DSMs and DTMs are prepared on the fly.

Perhaps we first create .vrt files for the DTMs and CSMs and pass to a custom ITS function?

My current workflow looks something like this

Iterate through tiles in parallel from 1 to N (parallel package)

read in lidar with readLAS within a parallel node, also read in dtm and dsm using raster(...)
use lastrees() with the CSM and watershed() to assign trees to points
use lasfilter to remove points not associated with trees
normalize tree heights and remove points below 2m with lasnormalize() and lasfilter()
generate hulls using tree_hulls() and compute .stdmetrics
replace values in PointSourceID with treeID*
write lidar to new .laz file with only tree data
write hulls to .shp file
write hulls attributes to .csv file (slightly redundant)

Perhaps a simple improvement without using LAScatalogs would be to extract data from a catalog file instead of reading in las files directly. This would enable me to buffer tiles slightly and remove edge artifacts (additional workflow steps involved), but would likely add some additional processing time.

My current thought is to build a new function and supply it to catalog_apply(), something like what is below, although in version 3.0 of lidR apparently the segment_trees() function can also accept a catalog.

tree_fn <- function(
 las_chunk
 , vrt_dsm
 , vrt_dtm
 , th
 , fn_metrics = .stdmetrics
 ,dir_las_out
 , dir_ply_out
 ,dir_csv_out
){
  #something like this:
  dtm_all  <- raster::raster(vrt_dtm)
  csm_all  <- raster::raster(vrt_dsm)
  ht_chunk   <- lidR::normalize_height(las_chunk, dtm_all)
  algo_all <- lidR::watershed(csm_all, th = th)
  ht_ws_chunk  <- lidR::lastrees(ht_chunk, algo_all)
  trs_chunk <- lidR::lasfilter(ht_ws_chunk, !is.na(treeID))
  hulls_chunk  <- lidR::tree_hulls(trs_chunk, type = "concave", concavity = 2, func = fn_metrics)

  #these steps I don't know how to do, or perhaps return and let catalog_apply handle writing somehow?
  if(F) lidR::writeLAS(trs_chunk,file.path(dir_las_out,"??"))
  if(F) rgdall::writeOGR(hulls_chunk,file.path(dir_ply_out,"??"))
  if(F) write.csv(hulls_chunk@data,file.path(dir_csv_out,"??"))

  return(hulls_chunk)
  #or return both?
  #return(list(hulls_chunk,trs_chunk))
}


texts: Place trace flags in bulk on selected features

texts:
I'm working on our company's geometric network looking for disconnected features. I have 102 purchase stations that I'm needing to trace from and was wondering if there was a quick way to select all and place flags on them?

I know I could run the Trace Geometric Network tool from toolbox, but I'd rather use the utility network analyst toolbar since the drawing method works quicker than a selection.

I've found some old threads where I'd have to write some code, but they weren't in python, and I'm not a strong programmer.


texts: Drawing arc by centerpoint, radius and azimuth in KML?

texts:
How do I generate lon/lat coordinates for an arc/sector?

Coordinates can be used/filled into a KML file. Pointers to geometric formula or a program applying the formula are sought.

Radius is 2 km and location is Europe North East. Accuracy is not a issue at this point.


texts:
I use QGIS 2.14.10

I am trying to clip a raster with a vector layer mask (shapefile)using the clipper tool (gdal). There is an offset betwween the resulting raster and the vector mask (see picture below).

I also tried with the saga script : clip raster with polygon. This one is working correctly but can only open a temporary file. If i try to save hard file i have this message : "ERROR 1: Can't load requested DLL: C:\PROGRA~1\QGIS2~1.14\apps\saga\dll\gdal_MrSID.dll
127: La procédure spécifiée est introuvable."

I need to clip a raster with a polygon ina script.


texts:
When opening an edit session Arc is telling me that some of my shape files' coordinate systems do not match the data frame.

Here is the data frame coordinate system:

  Stereo_70

  Projection: Double_Stereographic

  False_Easting: 500000.000000

  False_Northing: 500000.000000

  Central_Meridian: 25.000000

  Scale_Factor: 0.999750

  Latitude_Of_Origin: 46.000000

  Linear Unit: Meter

  GCS_Dealul_Piscului_1970

  Datum: D_Dealul_Piscului_1970

And here is the coordinate system that all the shape files in question have:

  Projected Coordinate System:  Stereo_70

  Projection:   Double_Stereographic

  False_Easting:    500000.00000000

  False_Northing:   500000.00000000

  Longitude_Of_Origin:  25.00000000

  Scale_Factor: 0.99975000

  Latitude_Of_Origin:   46.00000000

  Linear Unit:  Meter

  Geographic Coordinate System:  GCS_Dealul_Piscului_1970

  Datum:    D_Dealul_Piscului_1970

  Prime Meridian:   Greenwich

  Angular Unit:     Degree

Now is it just me or is anyone else having difficulty seeing the difference?  Doe this mean I'm going to have to reproject the files in question FROM their original projection, TO their original projection?  I'm not even sure that's possible but will try now.

And short of that I guess I just wipe the projections and define them again.


texts: Getting latitude and longitude of projected point using ArcPy?

texts:
I have a point feature in a feature class that is being accessed by ArcPy.  The point is projected but I need to find an efficient means to get the unprojected latitude and longitude for this point.

Is there a method other than reprojecting (unprojecting), getting a search cursor on the new feature class, finding the feature, then getting the lat/lon off the feature's shape?


texts: How to georeference a BC Land Survey in ArcGIS Pro

texts:
I am relatively new to ArcGIS Pro. I have a map with land parcels (from ParcelMap BC; let's call it the PMBC layer). I have added a DWG file from a BC land surveyor (let call it the BCLS layer). The land parcel boundaries do not align with each other.
Attempt 1: I tried applying control points under the CAD Layer Georeference tab to move the BCLS layer to align with the PMBC layer. The BCLS DWG file moved, but did not align with the PMBC layer.
Attempt 2: I looked at the PMBC layer properties and noticed that the PMBC land parcel layer was in 'NAD 1983 CSRS BC Environment Albers' so I used the Project tool to convert it to 'NAD 1983 UTM Zone 10N' to match the map and the rest of the layers on the map, including the BCLS layer. This resulted in a new protected layer located in the exact same position as the old layer.
Any tips on how to align the BCLS layer with the PMBC layer?