texts: Where did the texture fill option go in QGIS 1.9 Master?

texts:
In the previous version (1.8) of QGIS you could select to fill a polygon using a PNG or JPG image.  This option seems to have gone from QGIS 1.9 Master.  I can get it by selecting the "Old Symbology" button on the Style tab.  That's OK but surely this option has not been deprecated (or has it)?

Can anybody tell me where the fill by image option has disappeared to because the only other option I appear to have is SVG fill (SVG is very nice, but I also want the flexibility of using a PNG or JPG is that is what I have to hand) - or perhaps the option only exists as the 'old symbology'.  Also, does anybody know if it is intended to only allow filling polygons with SVG as opposed to PNG/JPG in the future?

Thanks


texts:
I have been using QGIS to insert point data accompanied by information for each point into PostGIS. When I open PGAdmin to take a look at the table rows I noticed a column near the start of each row called "geom" which is obviously point geometry with I believe a SRID of 4326.

Now my main questions are:

Does QGIS add this to the database or is this something that PostGIS adds automatically.
I'm going to be using a form which calls upon PHP to insert values into this databse table. Do I need to figure out the geometry or will PostGIS handle this for me?

I have looked at alot of information about SRID but a good bit of it is over my head.


texts:
I have a point where the objcet is and another point for mylocation.

I want to calculate the distance

>> object_location
<Point object at 0x7f9ceaac7be0>

>> mylocation
<Point object at 0x7f9cf9311b50>

>> mylocation.distance(location)
0.13980702957351368

But the result is wrong is there some different way with geodjango to calculate the distance between two points

This way in ipdb is correct

ipdb> Item.objects.filter(id=self.object.id).distance(mylocation, field_name="address__location")[0].distance
Distance(m=15368.0088756)


texts: How to keep the same style of a OSM map converted to SHP file

texts:
I saved a OSM file as a Shapefile, this is fixing my previous question how to handle bad layers.

My question is how to retain the original style from OSM plugin to SHP files, when I copy the style and try to paste it it says "Unknown renderer".

I mean, is it possible keep the same look I have from OpenStreetMap in the new SHP file?


texts: Emulating Polygon Neighbors tool in ArcGIS Desktop 10.0?

texts:
I'm using ArcGIS 10.0 SP5.

Have a set of 4000 polygons (featureclass A) and I'd like to summarise the values of the polygons (featureclass B) which touch the boundaries of these polygons (ie. those which neighbour A). The Summary statistics should include values in A (the field that I'd like to summarise is common to both featureclasses).

There is a tool in 10.1 and above called 'Polygon Neighbors' which would do exactly what I need, but that's not available in 10.0: https://pro.arcgis.com/en/pro-app/tool-reference/analysis/polygon-neighbors.htm

Is there a way to do this in 10.0?

Dissolving B polygons and using summarising statistics within this tool  will not work because some of the B polygons share a boundary with more than one A polygon. I need a separate summary for each A polygon.

This is not a duplicate question. The question which was marked as 'the same' as this asks where the 'Polygon Neighbors' tool is in 10.0. I realise that this tool isn't available in 10.0, I am asking how to produce the same output as this tool.


texts: Change the size of an OpenLayers.Control.OverviewMap

texts:
I create my OpenLayers.Control.OverviewMap with a given size but I would like to change it when certain conditions are met. I have tried two approaches with no luck...

My first try was just changing its width and height via CSS attributes:

$("div.olControlOverviewMapElement>div.olMap").css('width',320);

It does change its visual size but then the overview map doesn't work properly, the "new gained" space in the overviewmap is not usable (can't be used to drag the view rectangle)

My second try was simply changing its size property:

ovMap.size = new OpenLayers.Size(320,240);

but it doesn't have any effect to the overview map.

Is it possible to change its size once it's been created?

Extra information: Using OpenLayers-2.10, and in both cases I have tried calling the update() method of the control.


texts: PostGIS doesn't automatically convert WKB to GeoJSON

texts:
I discovered that somewhere in the upgrades from PostgreSQL 9.5 and PostGIS 2.2 to PostgreSQL 12.1 and PostGIS 3.0 the behavior for transforming a row  into JSON changed.

I am using the following example command:

SELECT row_to_json(tmp) FROM (SELECT '0101000020E61000006666666666665F403333333333332440'::geometry) tmp;

In 9.5/2.2  it produces the output:

                            row_to_json
-------------------------------------------------------------------
 {"geometry":"0101000020E61000006666666666665F403333333333332440"}
(1 row)

But 12.1/3.0 produces the output:

                       row_to_json
----------------------------------------------------------
 {"geometry":{"type":"Point","coordinates":[125.6,10.1]}}
(1 row)

Is there a setting to retain the functionality of 9.5/2.2 and not automatically convert my WKB to GeoJSON?

Edit:

The following command produces the error cannot cast type geometry to json in 9.5/2.2 and the expected GeoJSON value in 12.1/3.0:

SELECT '0101000020E61000006666666666665F403333333333332440'::geometry::json;

I am guessing it is related to how PostGIS now knows how to convert between the types, I wonder if there's a way to force it to interpret everything as text.


texts: QGIS printcomposer add reference to preset

texts:
When I make a visibility preset, I can use it in the template composer to 'lock' the layer styles.

I was wondering if I could use it as a variable in a textbox to display the name of the preset?


texts: How to find the minimum and maximum value of a region and give these values to the visualisation parameters in Google Earth Engine?

texts:
What I mean is: We are giving the min and max value to the visualization parameters but in my case min and max value changes according to geometry so how can i find the min and max value and give these values to the visualization parameters.

My code is:

var median = coll.reduce(ee.Reducer.median()).clip(geometry);

var extrema = median.select('CO_column_number_density_median').reduceRegion(ee.Reducer.minMax(),geometry,250);

var min1 = extrema.getNumber('CO_column_number_density_median_min');

var max1 = extrema.getNumber('CO_column_number_density_median_max');

  var visParams={
    min:min1,
    max:max1,
    opacity: 1.0,
    palette: ["blue","green","red","cyan","black","yellow","purple"],
    bands:['CO_column_number_density_median']
    };

  Map.addLayer(median,visParams, 'true-color composite');


texts: Cartographic techniques for symbolizing routing data

texts:
Some common problems with displaying routes created from street data, such as bus routes, include:

Routes may overlap themselves, such as when a route takes the same road in opposite directions
Routes may overlap other routes that take the same street.
Routes are often overly detailed, or "noisy", owing to their derivation from street data.
Indicating a route's direction unambiguously without overloading it with arrows can be a challenge.

What are some techniques for dealing with these issues in order to create a practical and usable route map, such as this Manhattan bus routes map? Do you know of any such public transport maps that were produced by automatic techniques rather than manual illustration/digitization?

Specifically I am using ArcGIS 10, but solutions involving other software are welcome as well.

One technique I have had some success with for dealing with 1) is using ArcMap's cartographic line symbols with a negative offset -- assuming the lines are oriented in the direction of travel, this will offset the lines such that the right and left-hand sides are distinct and oriented as you would expect (at least for right-hand traffic countries, for left-hand traffic countries one would use a positive offset). However the result tends to be a bit wonky around turns, interchanges and and U-turns, and is inflexible when dealing with multiple routes.

Another possibility that I am not very familiar with but which seems promising is cartographic representations -- however, I have not had much luck in finding examples of using representations for transportation maps and whether it might fit the bill.


texts: PostGIS ERROR: BOOM! Could not generate outside point!

texts:
This may be a continuation of a closed thread.

I want to sum the areas of polygons over categories in a large data set (avian species ranges), like this:

select sciname,
       sum(ST_Area(the_geom::geography)) / 1000000 as km2
from gis.aviandists_pl
group by sciname
order by sciname;

Some rows trigger 'ERROR: BOOM! Could not generate outside point!'. I've isolated a row that is causing this problem, but haven't found a way to upload it here. Based on the thread linked to above, I thought it might be a globe-spanning polygon, but it doesn't appear that way (see image below). This represents a single row from a geometry table, which surprises me because I would have thought that the multiple polygons that I'm observing would have to be represented by multiple records. Could this actually be a global polygon with small holes? More importantly, how can I work around this? The query covers over 10k rows, and I have no idea how many may pose this problem or how to find them. I'd be happy to share a shapefile of the problem data if that would help, but I haven't found a way to attach one here.


texts: QGIS: graphical-modeler to calculate area of polygons

texts: I'm trying to calculate the area of polygons using the field-calculator in QGIS (2.0.1). If I do this manualy it works fine, but if I do this using the graphical-modeler I get NULL as area for every Polygon. I just typet $area in the Formula field. Field-type is integer and length is 10.


texts: Need gps trackers for social science research

texts: I'm interested in studying people's movement over time (24 hours), and would like to do so with devices that also measure elevation very accurately. I understand there's stuff out the where I can buy these for about $20 a piece, as I need about 30 of them. The data will be used in arcgis. Can someone recommend a specific brand to buy?


texts: Quickly make server side spatial filters for WFS layers based on other layer's selected geometries in QGIS

texts:
In QGIS 2.12.1, I have a WFS layer which is huge so I don't want my users to load it entirely but only parts of it based on intersections with geometries.

Today, I'm using the query builder and execute filters like :

intersects(  $geometry , geom_from_wkt( 'Polygon ((29085.2677694238955155 5701652.79282033536583185, 23973.40924414427718148 5683250.10212932899594307, 39820.17067251115804538 5683250.10212932899594307, 37264.24140987131977454 5698074.49185263924300671, 37264.24140987131977454 5698074.49185263924300671, 29085.2677694238955155 5701652.79282033536583185))'))

I need to use the plugin getWKT and copy/paste the WKT which is not practical for beginners. Also this plugin doesn't handle multiple selected features.

So I'd like to know if it is possible to easily do the same job with geometries from selected features of another layer?

EDIT :
Inspired by this article I've created a user defined function :

from qgis.utils import qgsfunction, iface
from qgis.core import QGis
@qgsfunction(0, "Python")
def getSelectedFeaturesInWKT(values, feature, parent):
  """
    Returns the wkt of selected features
  """
  layer = iface.activeLayer()
  features = layer.selectedFeatures()
  geom = features[0].geometry()
  return geom.exportToWkt()

This function return the first (to test then I want to get all features) selected feature in wkt.

Then I load this script and modify my query :

intersects(  $geometry , geom_from_wkt(getSelectedFeaturesInWKT()))

The expression is valid but the layer created has the URI of the WFS server without the GML filter.

I have tested getSelectedFeaturesInWKT() alone and it works.


texts: Problem with previewing layer on geoserver ( instead of downloading wms file )

texts:
When I hit OpenLayers to preview my layer in Geoserver, instead of opening a new tab with OpenLayers it automatically downloads a wms file
I just opened the wms file and says :

<?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE ServiceExceptionReport SYSTEM "http://localhost:8080/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd"> <ServiceExceptionReport version="1.1.1" >   <ServiceException>
  java.lang.NullPointerException


texts:
I have two rasters with different extents that overlap partly:

ext1 <- extent(99500,700500,249500,600500)
ext2 <- extent(1356,643990,269903,690810)

#create rasters
rr1 <- raster(matrix(runif(1000000,0,100), 1000, 1000))
rr2 <- raster(matrix(runif(1000000,0,100), 1000, 1000))

extent(rr1) <- ext1
extent(rr2) <- ext2

How can I plot both rasters in one plot?

I've tried

ext3 <- extent(min(ext1[1],ext2[1]),max(ext1[2],ext2[2]),
               min(ext1[3],ext2[3]),max(ext1[4],ext2[4]))
plot(rr1, ext = ext3)

But the only the extent of rr1 is used.


texts:
I recently updated R and all of it's packages, and my code for reading in the layers of a geodatabase no longer works.
It works fine on my windows desktop, but doesn't work on my macbook (I change the filepath for the appropriate machine).
I get the same behaviour with any shapefile or geodatabase at any file path.
I'm essentially following the instructions listed here: https://stackoverflow.com/questions/52340442/cannot-open-data-source-gdb-in-r
r1_gdb <-  "/Volumes/LaCie/2019 PRISM GIS clean/Region1_Cleaned.gdb"

subset(ogrDrivers(), grepl("GDB", name))
          name    long_name write  copy isVector
50 OpenFileGDB ESRI FileGDB FALSE FALSE     TRUE

r1_lyrlist <- ogrListLayers(r1_gdb)

Error in ogrListLayers(r1_gdb) : Cannot open data source

What strikes me as weird is that the results of other functions leads me to believe this should be working fine.
file.exists(r1_gdb)
[1] TRUE

list.files(r1_gdb)
 [1] "a00000001.freelist"
 [2] "a00000001.gdbindexes"
 [3] "a00000001.gdbtable"
 [4] "a00000001.gdbtablx"
.....

Here is info about my version of GDAL on my macbook.
> library(rgdal)
rgdal: version: 1.5-19, (SVN revision 1092)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.1.4, released 2020/10/20
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/proj
Linking to sp version:1.4-4
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
> getGDALVersionInfo()
[1] "GDAL 3.1.4, released 2020/10/20"

It works fine on my windows desktop. It seems I have a different version of GDAL on the two computers. But it's the windows computer with the older version of GDAL where the code doesn't work.
> library(rgdal)
Loading required package: sp
rgdal: version: 1.5-21, (SVN revision 1105)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
Path to GDAL shared files: C:/Users/andersonch/Documents/R/win-library/4.0/rgdal/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
Path to PROJ shared files: C:/Users/andersonch/Documents/R/win-library/4.0/rgdal/proj
Linking to sp version:1.4-5
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

Any suggestions as to what I should investigate?


texts:
How is it possible to perform multiple processing steps for one layer without creating additional ones?

I wanted to add several fields to attribute table using qgis:fieldcalculator but every time it works new layer is created. Is it possible to redirect output from first .runalg to second and so on until getting final result without additional layers?

first = processing.runalg('qgis:fieldcalculator', 'S1_PROJ_1992[PTS]', 'X_2000', 0, 10, 1, True, '@row_number', None)
second = processing.runalg('qgis:fieldcalculator', first['OUTPUT'], 'Y_2000', 0, 10, 1, True, '@row_number', None)
processing.runandload('qgis:fieldcalculator', second['OUTPUT'], 'Y_2000', 0, 10, 1, True, '@row_number*(-1)', )


texts: How to view attribute table in layer within ArcGIS Server map service?

texts: I'm new to ArcGIS Server and am trying to publish layers within an mxd to ArcGIS Server. When I bring the published map service back into Arc there is no Open Attribute Table option although I am able to identify the data. Is there away of displaying the attribute table?


texts: Why does QGIS's Clipper tell me "Computed -srcwin falls outside raster size" and how do I address it?

texts:
I'm trying to clip a heightmap to a small area that I'm interested in. When I attempt to use Raster > Extraction > Clipper I get an error dialog with the message ""Computed -srcwin falls outside raster size of..." with the dimensions.

I tried to select a smaller area but continue to get the error message and am not sure what it's complaining about and how to resolve the issue so I can clip the heightmap.

What am I doing wrong or need to adjust?


texts: Using Craig Retro-azimuthal (Mecca) projection in QGIS?

texts:
I'm using QGIS 3.4.12 and wanted to make a map with the Craig Retro-azimuthal projection.

The projection is not implemented in Proj4 and I couldn't find information about implementing new projections in the QGIS documentation. Is there anyway to "Install" the projection in QGIS?


texts: Dissolve in QGIS to merge several processes in different regions into one process with all values

texts:
I am working on a dataset where I combined geographical information of cities in Europe and changes of landcover in these regions. I am working with the CORINE Dataset and I would like to know about the variation of changes in the whole area in different time periods. In the CORINE dataset you can find a lot of different changes (e.g. dense urban fabric -> parc areas (112-212) etc.). I'd like to combine all changes which are similar (see example above) in the observed area.

I managed to do a dissove where the output was only one field in the attribute table for one change (e.g 112-212) and I thought that's it. After some troubles I controlled the results and saw that the information about the size of the area is not the size of all combined areas what makes it unusable for me.

Some pictures might explain more than my words:

2 Datasets:

I merged the two datasets to reduce information and to have the changes just at the area where the cities are and then used Dissolve to field "CHANGE" (worked) but I also need the cumulated information about all the areas where that change happened (area h and the number of areas which are now packed into one line in the attribute table (did not work)).

Do you have any solution for my problem?


texts: Calculating total PV potential of roof using QGIS

texts:
I want to calculate total PV potential of a roofs in a certain area. To do that I made a dsm, ran SAGA "Potential Incoming Solar Radiation" algorithm (just playing with it only for now - I need to work out how to obtain proper values) and cut the output raster to vector with a roofs outlines.
As a result I received nicely cut raster that fits to size of each roof, as visible on image below.

What tool should I use to obtain exact value of irradiance for area of individual roofs?
To be more accurate: not only visualize distribution of irradiance on the roof but exact number of kWh/sq m for given area?


texts: Calculating spacing of polylines along transects using ArcPy?

texts:
I have a shapefile with a number of roughly parallel lines in it, and I want to calculate the spacing between these lines using a number of transects across the map.

How should I go about this? I have a bit of experience writing ArcPy scripts, so would be happy to write some code to do this, I just can't quite work out where to start. My best guess would be to draw the transect as a line, somehow create points where the transect line intersects with any other lines, and then calculate the average distance between the points. Does this seem like a sensible approach? If not, what would you suggest?

I need to be able to do this scriptably (from a Python script using ArcPy), which means I won't have access to functions that are only available from toolbar buttons.


texts:
I am using the marine boundaries data from MarineRegions.org and I'm trying to use cylindrical equal area with a -160 central meridian.  I successfully removed the dateline from polygons basically using the method outlined in this GISSE question.  However, the marine boundaries dataset the polygons are singlepart, even with the split and ArcGIS is not allowing me to explode to multi-part (to then dissolve)... the option is grey-ed out.  Can anyone provide a solution to the polygon-crossing-international-dateline problem when the original polygon is singlepart?

Basically, is there a way to accomplish this in ArcGIS?


texts: Does "size scale field" in QGIS actually work?

texts:
I'm having mixed results with the 'size scale field' aka graduated symbols in QGIS.

Some forum posts out there say this feature is broken - others can get it to work.


texts: Field Join by Approximate Data in ArcMap?

texts:
I've run into a bit of an issue. I recently pulled data from the CDC website concerning Suicide rates by county. However, the csv file does not have locational data. Luckily, my project is about the correlation between suicide rate and % of abandoned homes in the nation (by county).

The abandoned homes shapefile does have locational data and I wanted to do a table join based on county name or county code; however, for both instances, the data within the field is just a little different.

For instance, Autauga County, AL versus Autauga County, Alabama.
As you can see, the names almost align. The default method of joining the fields based on county names brings up a lot of "Null" values.


texts: Invalid URL error when using suggestions in ESRI locator service

texts:
I have created and published a locator service with suggestions turned on. I enabled suggestions when building the locator files and after publishing them as a service. The service published correctly and geocoding from the rest endpoint works fine. E.G.

http://localhost:6080/arcgis/rest/services/BLDGADDRESS/GeocodeServer/findAddressCandidates?Street=&Single+Line+Input=88+quale&category=&outFields=*&f=html

this works fine

However using the suggest option give an invalid URL error e.g.

http://localhost:6080/arcgis/rest/services/BLDGADDRESS/GeocodeServer/suggest?text=555+night&f=html

gives a 404 bad URL error. Any ideas why?


texts: Adapt attribute table QGIS legend from bottom

texts: I have an attribute table controlled by Atlas on map composer from QGIS. I need to adapt the table for each map and its features but from bottom. In other words, I need to make grow the table from bottom to top. By default table grows from top to bottom


texts: Creating polygon based on spatial extent in ArcGIS Desktop?

texts:
Is there anyway that I can input:

The Coordinates System (projected is OK)
The bounding box

and create a rectangle polygon based on these?


texts:
I am new to mapping in R.

I want to have a density map of US for number of phone calls for different areas. Each area includes some zip codes. I only have number of calls for each zip code and I want to show it on a map of US by using zip codes. As I don't have polygon coordinates for my areas, I decided to use the zip3 shapefile which is available on the web and merge my data with the existing shapefile.

To be more clear, my area is formed by zips 006, 010 and 011. I defined an ID for each area and if different zips are related to the same area, they will have identical ids. Hence, my Excel file includes area name, zip3 and ids related to each zip3. I used zip3 from my Excel file as a joining variable to the shapefile (shapefile has a zip3 column) and allocated different ids to zip3s which are in the shapefile.

At this time, I wanted to merge polygons of zip3 s which forms an area to build a map of my areas from the shapefile of zip3s. I believe this is the right approach but it is not working. My problem is that maybe I'm reading my shapefile to a wrong data frame. I used readShapePoly() to read my shapefile.

My code is as follows:

data <- read.csv(“zipmatch.csv”,colClasses=c(“character”,”character”,”numeric”),header=TRUE)
data2 <- readShapePoly(“zip3.shp”)
data3 <- merge(data2,data,by="ZIP3")
final <- unionSpatialPolygons(data3,"id")

It gives me the following error at this time:

  Error in unionSpatialPolygons(data3, "id") : input lengths differ


texts: Create feature class with ArcObjects in SDE geodatabase returns No editable Layer

texts:
I Want To create Feature Class in an enterprise geodatabase with ArcObjects. This is the small part of the our plugin. In our case, we want to Create Feature Class in the database Identical to the layer that is in TOC. It creates that successfully, but when we add it to the TOC and want to edit it, ArcMap Return the message that "This data is not editable". The user has full permissions .
It is possible to create Feature Class manually in geodatabase With database Connection in Catalog tree; It successfully created and it is editable and i can Register it as versioned . But When i try to create feature Class with Arcobjects it returns No editable Layer and i cant edit feature Class.
  SQL Server in use: Microsoft SQL Server 2008 Express Edition (64-bit) ,  version of  ArcGIS in use: 10.2.2, the license level of ArcGIS: Advanced (ArcInfo) Concurrent Use.

This is our code:

   Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.SqlWorkspaceFactory");
        IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance
            (factoryType);
        IPropertySet connectionProps = new PropertySetClass();
         connectionProps.SetProperty("SERVER", "gSERVER");
        connectionProps.SetProperty("INSTANCE", "sde:sqlserver:gSERVER-pc");
        connectionProps.SetProperty("DATABASE", "Ga");
        connectionProps.SetProperty("DBCLIENT", "SQLServer");
        connectionProps.SetProperty("SERVERINSTANCE", "gSERVER-pc\\sqlraSERVER");
        connectionProps.SetProperty("AUTHENTICATION_MODE", "OSA");
        IFeatureWorkspace workspace = (IFeatureWorkspace)workspaceFactory.Open(connectionProps, 0);
        IMxDocument mxd = ArcMap.Application.Document as IMxDocument;
        IMap map = mxd.FocusMap;
        ILayer layer = map.Layer[0] as ILayer;
        IFeatureLayer featurlayer = layer as IFeatureLayer;
        IFeatureClass fc = featurlayer.FeatureClass;
        IWorkspace wrk = workspace as IWorkspace;
        IFields fields = featurlayer.FeatureClass.Fields;
        IFeatureWorkspace hhj = workspace as IFeatureWorkspace;

        IFeatureClass outputFeatClass = hhj.CreateFeatureClass(featurlayer.Name, fields, featurlayer.FeatureClass.CLSID, null, featurlayer.FeatureClass.FeatureType, "SHAPE", "");


texts: If the GPS navigation message takes 12 1/2 minutes to cycle, how can receivers update your position every second?

texts:
According to Wikipedia's article on GPS,

  Each GPS satellite continuously broadcasts a navigation message... Each complete message takes 750 seconds (12 1/2 minutes) to complete.

My understanding of how GPS works is that the receiver receives the broadcast message and calculates location based on that. So how can a GPS system start delivering reasonably accurate locations within a few seconds?

Wikipedia mentions that the message is divided into frames and sub-frames:

  The message structure has a basic format of a 1500-bit-long frame made up of five subframes, each subframe being 300 bits (6 seconds) long... [A] complete data message requires the transmission of 25 full frames.... [T]his gives 750 seconds to transmit an entire almanac message (GPS). Each 30-second frame begins precisely on the minute or half-minute as indicated by the atomic clock on each satellite.

But the iPhone's GPS, for example, delivers about a location point a second. Even if a frame or a subframe is enough, that's 6-30 seconds. How can the hardware update its position with reasonable (reported) accuracy every single second? Is the iPhone lying to me?


texts: Summing multiple image outputs in earth engine?

texts:
I am doing some GDD data analysis for a project and am having trouble coming up with a way to effectively analyze a specific time of the season over multiple years.

What I am trying to do is find what areas accumulate GDD and what areas don't during a critical period of fruit production over a period of multiple years.

I have the code below (which does the calculations for GDD), but i am stuck on how to run that calculation for each filter date that I have set (04-15-xx, 05-15-xx) and then either sum or average the image outputs.

This works for an individual year (by replacing ".map(ClipToCol)" with the filterDate function, but i want to average or sum multiple years.

Ideas?

// This function clips images to the ROI feature collection
var clipToCol = function(image){
  return image.clip(geometry);
};

var dataset = ee.ImageCollection('MODIS/006/MYD11A1')
                               .map(clipToCol);

var landSurfaceTemperature = dataset.select(['LST_Day_1km', 'LST_Night_1km']);

// According to: https://gis.stackexchange.com/questions/307548/getting-temperature-data-of-given-point-using-modis-lst-data
// map over the image collection and use server side functions
var tempToFahr = landSurfaceTemperature.map(function(image){
  var props = image.toDictionary(image.propertyNames());
  var Fahr = (image.multiply(0.02).subtract(273.15)).multiply(1.8).add(32);
  // Mask where one of Night or day temperature has no data
  var FahrMasks = Fahr.updateMask(Fahr.select('LST_Day_1km')).updateMask(Fahr.select('LST_Night_1km'));
  // we assume that the night temperature is the min temp, and the day temperature is the max temperature
  var meanFahr = FahrMasks.reduce('mean').rename('meanTemp');
  // Calculate the GGD and make all the negative values 0 (see: https://en.wikipedia.org/wiki/Growing_degree-day)
  var GDD = meanFahr.subtract(52.5).rename('GDD');
  var GDDnonNeg = GDD.where(GDD.lt(0), 0).rename('GDDnonNeg');
  return ee.Image(Fahr.addBands([meanFahr, GDD, GDDnonNeg]).setMulti(props));
});

// calculate the sum of GGD values
var summed = tempToFahr.select('GDDnonNeg').sum().rename('summedGDD');

var landSurfaceTemperatureVis = {
  min: 0,
  max: 100,
  bands: ['LST_Day_1km'],
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};


texts: Remove last line segment of an OpenLayers.Control.Measure

texts:
i'm using OpenLayers library, i'm trying to delete the last line segment of a drawing realized with the following code:

layerLine = new OpenLayers.Layer.Vector("layerLine");
layerLine.id = "id_layer_line";
map.addLayer(layerLine);

lineControl = new OpenLayers.Control.Measure(OpenLayers.Handler.Path, {
        id: 'id_layer_line',
        persist: true,
        layer: layerLine ,
        //immediate:true,
        handlerOptions: {
            layerOptions: {styleMap: styleMap},
        }
    }
);

But the layerLine has no features (empty array), it means i can't access to line segments, is there are any way to do this?


texts: Establishing DGPS Connection With Old Equipment?

texts:
I am trying to do a differential surveying project with some outdated equipment, specifically:

-An old AllegroCX personal handheld computer (Win 4.1) with SurvCE loaded on it to collect data and establish communication

-Two Novatel ProPak V3 Recievers (one for base, one for rover)

-Two external radio communicators

-Two pinwheel style antennas

I am able to establish a connection with the rover's Novatel reviever and AllegroCX using a serial cable, but I cannot get it to find the base receiver using radio communications. Is this possible wirelessly, or do I need to set up the base station via serial cable as well? Do I need to have a PC connected to the base receiver?


texts:
I have a water tank data(.shp) and I need to label it using field NAME. However I just want to label water tank that have names like "Tangki Air Sembulan". If in the field NAME only states "Tangki Air", I do not want the label to appear. Im using Arcgis 10.2 and label using Maplex. I try to use SQL Query like below:

NAM NOT Like '%Tangki Air%'
AND NAM NOT Like '%TANGKI AIR%'

but it's not working. I cannot edit the field name so is there a way that can help me solve this problem?


texts: Openlayers: How to highlight a point on mouse over a line

texts:
Any suggestions on how to highlight a point on mouse over a line, I've read points from a javascript array, lines are drawn by using them and I want to highlight some points in that array by moving the mouse pointer through the line and near to that point. However, i'm  stuck on actually how to get started with this.

Any guidance and help is appreciated..!


texts: changing attribute table automatically in map composer

texts:
I would like to add an attribute Table to the map composer, which so far is working, I also managed to display it as two coloums. I'm using QGIS 2.18.3 (Las Palmas) for a project where I have about 15-20 maps with the same extent, each displaying a different attribute though. I am able to use the atlas to create a seperate map for each of the 15-20 attributes (through visibility Presets), but I am not able to automatically change the attribute table accordingly in my map composer.
The maps composer switches through the different maps, but always displays the attribute table of the first map. I'm not sure where I'm going wrong.


texts: Struggling with in_memory to erase faster

texts:
I want to delete all points inside a polygon layer:

import arcpy
import time

arcpy.env.workspace = r"E:\DensityMaps\DensityMapsTest1.gdb"
arcpy.env.overwriteOutput = True

start_time = time.time()

arcpy.Erase_analysis("Test_Points_Project", "Europe_coastline_NoLakesRivers", "PointsPolyErase1")

print "Erase ended "+("--- %s seconds ---" % (time.time() - start_time))

This lasts 36,9 seconds

Then I read this post saying that the process can be dramatically optimized using in_memory. I give it a try:

import arcpy
import time

arcpy.env.overwriteOutput = True
start_time = time.time()

#Copy points to in_memory
inFeature = r"E:\DensityMaps\DensityMapsTest1.gdb\Test_Points_Project"
memoryFeature = "in_memory\PointsMemory"

arcpy.CopyFeatures_management(inFeature, memoryFeature)

#Erase points that are inside polygons
arcpy.Erase_analysis(memoryFeature, r"E:\DensityMaps\DensityMapsTest1.gdb\Europe_coastline_NoLakesRivers", "PointsMemory")

#Copy points back to directory
arcpy.CopyFeatures_management("PointsMemory", r"E:\DensityMaps\DensityMapsTest1.gdb\PointsPolyErase")

print "Erase ended "+("--- %s seconds ---" % (time.time() - start_time))

This lasts 38.4 seconds, 1,5 seconds slower than the previous script. I am obviously doing something wrong. Very wrong.

It's very important to make the process faster because the original point file contains 22 millions points (1,5 GB). If I don't manage to do it fast the operation can last about 24 hours.

Here are the files if you want to test


texts: List of methods of determining the projection of a shapefile when .prj is missing

texts:
Possible Duplicate:
  Best practices for identifying the unknown coordinate system of a Shapefile

It would be helpful to mention what information is needed to perform each technique. What are situations in which this may happen, and the best ways to deal with it?


texts:
In QGIS 3 there is the possibility to illustrate proportional data like population with circle of proportional size with the "assistant" option.

But once in the print composer section, there is an empty symbol in the legend.
How to show 2 or 3 circles of proportional size in the legend?


texts: Count number of overlapping map matched GPS trajectories with road network

texts:
I map matched a sequence of GPS-points to a road network using the Offline-Map Matching plugin for QGIS. This map matched trajectory passes certain roads segments multiple times (for example, the selected red matched trajectory segment in the picture below has two unique fid numbers, 60 and 359, see left).

In the picture below, only trajectory with fid 60 is selected. As you can see, these trajectory segments are generated between two consecutive GPS points.

In order to visualise where most vehicle activity occurred in the network (with the line thickness/size being proportional to the number of overlapping map matched trajectories), I need to analyse the number of overlapping map-matched trajectories per road link in the road network.

I have tried several things but unfortunately without success. Any ideas on how to do this in QGIS?


texts:
I need a polygon layer to draw in order based on a date field...  What I'm trying to do:

I need to display an easement feature class with two different symbols: one for easements and one for vacations.  In addition, I need this layer to draw in an order that represents the date in which they were recorded (stored in a date field).  Using symbol levels is not an option because I am working with thousands of different dates (Also, ArcMap symbology options do not recognize date fields as dates...).  The sort tool, which reorders the features in the database, is also not an option because this feature class is constantly being updated and it is not feasible for this tool to be ran frequently.

Any ideas?


texts:
I've got the following (2D!) situation of overlapping polygons:

I would like to have the polygon border dissolved such that I have only one polygon left. I thought this was a definitive case for the dissolve tool, but somehow the result is bizarre:

What did I do wrong? I haven't chosen the attribute to dissolve on so I thought then all polygons are dissolved and should bring the desired result.
Can anybody help?

EDIT:
A clipped version of the data can be found here:

https://www.dropbox.com/s/xm9ykgnrbleutgr/extracted_area.zip?dl=0


texts: My Loop Route in Google Transit is showing a 1 hour wait. Any help with my GTFS?

texts:
I have a Loop Route that's showing a 1-hour wait time at the last/first stop.

Google Query example

This is the last/first stop of my loop.

BLU,06:30:00,06:30:00,PRICE_ST,1,To Crescent Harbor
...
BLU,07:28:00,07:28:00,ARR_TRANSF,18,To Crescent Harbor
BLU,07:30:00,07:30:00,PRICE_ST,19,To Crescent Harbor

Have I misunderstood how to use the frequencies.txt? My understanding is that it adds 3600 seconds (1 hour) to the first stop's time and continues on from there.

trip_id,start_time,end_time,headway_secs,exact_times
BLU,06:30:00,19:30:00,3600,1

Or is there some issue with not enough time between the last stop and the first stop?


texts: I want to add another field to an attribute table, a number field. This number field will show in my opinion the experience visiting the ski resort, from 1 (the worst) to 10 (the best). Those that I have not visited will get 0 by default when I add the field, then I want to populate the ones I have visited with a value from 1 to 10. How do I go about this?


texts: Importing GDAL in spyder from a local GUI in cygwin

texts:
I have installed anaconda on linux (python 2.7).  Since gdal is not pre-loaded I load it with conda install gdal, which is fine.

When I do:

from osgeo import gdal

I get the following error though:

Traceback (most recent call last):

  File "<ipython-input-3-e4290214a986>", line 5, in <module>
    from osgeo import gdal

  File "/home/spotter/anaconda2/lib/python2.7/site-packages/osgeo/__init__.py", line 21, in <module>
    _gdal = swig_import_helper()

  File "/home/spotter/anaconda2/lib/python2.7/site-packages/osgeo/__init__.py", line 17, in swig_import_helper
    _mod = imp.load_module('_gdal', fp, pathname, description)

ImportError: libjpeg.so.8: cannot open shared object file: No such file or directory

I have done this same process on a windows with anaconda and gdal but no error.

My linux is Red Hat Enterprise Linux ComputeNode release 6.7 (Santiago)

EDIT:

Ok i think I am onto something.  What I am doing is using cygwin to enter a linux machine from a windows 7 computer. Within the linux machine if I do python -c "import gdal" all is fine.  The hang up is that I also want to use spyder linked to a gui in my local machine.  To get spyder to startup I need to change the PATH to export PATH=/home/usrname/anaconda2/bin:$PATH

as soon as I do this I must break some sort of connection with gdal.


texts: Calibrating routes gives interpolation error in ArcGIS Desktop?

texts:
I know that copy a question is not exactly great way to ask on this forum, but my situation is exactly the same as this unresolved issue on the (closed) ESRI forum:

  I have a stream centreline (fairly curvy) which can be converted to a
  route. I have a few points on that line which have stream bed
  elevation. The bed elevation varies along the stream, and does not
  necessarily increase as you go up the stream (e.g., goes up and down a
  little bit).

  I want to calibrate the route using the bed elevations, so that I can
  interpolate an elevation value between the known points.

  When I attempt to calibrate the stream route with the known points, it
  appears to work. But if I examine the measure values of individual
  vertices along the route, they are incorrect. Vertices that fall
  exactly where there is a known value (a calibration point) have the
  correct elevation, but vertices between the known points have
  incorrect values. The incorrect values are not 0, which makes me think
  it is calculating a value somehow. But it is not doing as I expect.

  For example, between known points of 3.11 and 3.56, you'd expect
  values between 3.11 and 3.56. Instead, I might get values of 3.78,
  3.78, 3.78, 3.79, etc.


texts: Is there an arcpy tool to delete/unregister a replica geodatabase?

texts:
Seems odd but I can't find an arcpy or geoprocessing tool to delete/unregister replicas.

When you create a replica it gets registered with the GDB that the data comes from. Afterwards you can't create another replica with the same name, you must first unregister the first one. To unregister a replica after it has been created you need to right click on the GDB in ArcCatalog, chose distributed databases and then manage replicas. From there you can pick an individual replica from all the replicas in that GDB and unregister it. This is what I would like to be able to do programmatically.

There are arcpy tools to create and sync replicas but not delete. Am I just being blind or does such a tool not exist?


texts: QGIS Shapefile - Units/Editing Shapefile

texts: I have a Shapefile that has properties/parcels as polygons. I am looking to select by parcels that are bigger than a certain size in acres. The problem is that the Shapefile has the properties/parcel sizes in BOTH square meters and acres. Some are in square meters and others are in acres, and there is a column that indicates whether the value is in acres (A) or if it is in square meters (S). Is there a formula that I can enter to create a new column that converts all of the numbers that are in square meters into acres?


texts: Are there any open source tools to detect and fill sinks on a DEM?

texts: Are there any open source or inexpensive tools to detect and fill sinks on a DEM? ArcGIS Spatial Analyst is just out of my price range.


texts: Multi Value User Input as Variable in SQL Statement

texts:
I'm attempting to write a script that takes a user's selection as a parameter and uses that as part of the SelectByAttributes function. The parameter represents field attributes and is multi value, and the user can select one or multiple.

So what I have at this point is the list of all available attributes, and each selected as the defualt.

Also, Here is the tool validation code and Parameter Properties:

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
    """Setup arcpy and the list of tool parameters."""
    self.params = arcpy.GetParameterInfo()

  def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""
    return

  def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    shape = "G:\_DPD\STEVEN\FIELD.shp"
    rows = arcpy.SearchCursor(shape)
    self.params[0].filter.list = sorted(list(set(r.getValue('DPD_FLD_TY') for r in rows)))
    all = self.params[0].filter.list
    self.params[0].value = all
    del rows
    return

  def updateMessages(self):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return

I want the user to have the option to run the tool with multiple attributes selected or
with just one selected, and use that as part of the SQL statement within the SelectByAttribute function. For example, if the user selects "GAS" and "UNKNOWN", the SelectByAttributes function would look like

arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",""""DPD_FLD_TY" = 'GAS' OR "DPD_FLD_TY" ='UNKNOWN'""")

So my main question is: what would the code be to check how many attributes are selected and based on that, have the SQL statement contain the proper amount of "OR" statements? Also, would the parameter be considered just one, or would it be multiple parameters i.e. GetParameterAsText(0) vs. GetParameterAsText(0)...GetParameterAsText(3)

I'm sorry if this is confusing, but I don't know exactly where to start with this problem.

This is how I am referencing the layer:

import arcpy

mxd = arcpy.mapping.MapDocument('CURRENT')
mm = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyr = arcpy.mapping.ListLayers(mxd, 'FIELD')[0]
field = arcpy.GetParameterAsText(0)


texts: Request Double Type Feature from Geoserver in JavaScript

texts:
I have set up an application according to this workshop (using postgis, pgrouting, geoserver, openlayers 3). In order to get a value for the total cost of each route selected I altered the SQL View in Geoserver in step 9 to:

SELECT ST_MakeLine(route.geom), sum(route.cost) AS totalcost
FROM (
    SELECT geom, cost
    FROM pgr_leastenergy_rec('ways', %x1%, %y1%, %x2%, %y2%)
    ORDER BY seq
        ) AS route

Everything worked out fine - I got this result in Geoserver:

Since the geometry is already implemented in the routing interface of the next step, I wanted to add the cost factor as well - as an output right under the map. I already tried a few things, unfortunately I am quiet new especially to Javascript. The viewparams variable and the attached request must be altered somehow. Any idea how to achieve that goal?


texts: Cloud cover percentage in Google Earth Engine for Sentinel-2

texts:
I'm trying to obtain the cloud cover percentage of my ROI in a specific time interval using Sentinel-2.

//Sentinel collection
var sentinel2 = ee.ImageCollection("COPERNICUS/S2")
.filterBounds(ROI)
.filterDate('2017-06-01','2018-06-01');

//Add a cloud cover band
var cloud = function(sentinel2) {
return sentinel2.addBands(sentinel2.metadata('CLOUD_COVERAGE_ASSESSMENT')); };
print(sentinel2.map(cloud));

var getCloudScores = function(sentinel2){
var value = ee.Image(sentinel2).get('CLOUD_COVERAGE_ASSESSMENT');
return ee.Feature(null, {'score': value})
};

var s2 = sentinel2.filterBounds(ROI);
var results = s2.map(getCloudScores);
var chart = ui.Chart.feature.byFeature(results)
    .setChartType('LineChart')
    .setOptions({
      title: 'Cloud cover',
      hAxis: {
        title: 'Date',},
      vAxis: {
        title: 'Clouds'
      },
      lineWidth: 1,
      pointSize: 0
});

print(chart);

So far this code has worked, but I was wondering if there is another way to get it?


texts: Event layer to feature class - fastest method?

texts:
There are three methods to convert the output of a Make XY Event Layer to a Feature Class:

Copy Features
Feature to Point
Feature Class to Feature Class

Which one is the fastest method? I need to convert a huge amount of points stored in a csv table (1,5 GB  file) and I wonder what's the best way to do it.

I will use Background geoprocessing on a 64-bit machine.


texts:
Working on ArcMap Desktop 10.8;
I have a point layer with a field called SITEROAD (string) from which I need to extract ONLY the first word before a space, and put it into a new field, STREETNAME.
Example:
Aberdeen way   needs to be calculated to Aberdeen
ACACIA WAY     needs to be calculated to  Acacia
Aegean WY      needs to be calculated to  Aegean

I am trying:
def getFirst(inField):

  if inField is None:

   return None

  else:

   return inField.Split('')[0]


texts: Python Geospatial development book error?

texts:
I am currently reading Python Geospatial Development by Erik Westra but I am stuck at the first hurdle where the book suggests that the easiest way to install gdal onto a windows machine is using FWtools from the website http://fwtools.maptools.org

The first problem is that this is a really old website.

It then says to check if it's successful type >>> import osgeo
this is where it fails and I get the error:

Traceback (most recent call last):
  File "<stdin>", line 1, in ?
ImportError: No module named gdal

The error seems to have come up a number of times on stack exchange:

Installing GDAL with Python on windows? (GIS.SE)
Unable to install Python and GDAL (DLL load failed) (SO)
can't import gdal in python? (SO)

with the answer being to use OSGEO4W.

The fwtools website hasn't been updated since about 2010 or 2011 and the book is from 2013 second edition so surely they had time to update the installation section, it seems to be a fundamental flaw in the book which has been rated highly on Amazon and by a number of people who know way more about python then I do, So my question is have I missed something really obvious in the book or that is considered general  python GIS knowledge? There is nothing in the errata on the books Publishers website Packt Publishing. I don't want to rate the book badly if it something I have missed.


texts: Adding 'freehand' data to raster mapping in QGIS

texts: Having just worked out how to import raster mapping into QGIS at the correct scale, I'd like to now add useful information to the map like marked routes, labels etc. How can I do this in a 'freehand' manner, rather than importing existing datasets?


texts: Finding most remote spot in Eastern United States?

texts:
I have this hunch that I could travel no more than 50 miles from the most remote spot in the Eastern United States (east of the Mississippi River), in the direction of the nearest road, and find a road.

Definitions: Most Remote: Spot furthest from a road. Road: Google Maps definition of a road.

How could I prove or disprove this claim  (i.e where is the most remote spot in the Eastern US)?


texts: 10 days composite from Sentinel 2 Images Collection

texts:
A monthly composite can be made using the following code, I take example from Sentinel 2 L1C:

var geometry = ee.Geometry.Point([109.061, 11.78]);
var col = ee.ImageCollection('COPERNICUS/S2')
          .filterBounds(geometry).filterDate('2018-01-01','2019-01-01');
//Monthly Images
var months = ee.List.sequence(1,12);
var years = ee.List.sequence(2018,2018);
// monthly composite
var monthlyImages =  ee.ImageCollection.fromImages(
  years.map(function (y) {
  return months.map(function(m){
  var w = col.select('B1').filter(ee.Filter.calendarRange(y, y, 'year'))
           .filter(ee.Filter.calendarRange(m, m, 'month'))
           .mean()
           //.median()
           //.min()
           .set('system:time_start',ee.Date.fromYMD(y,m,1));
  return w.set('year', y)
           .set('month', m)
           .set('date', ee.Date.fromYMD(y,m,1))
           .set('system:time_start',ee.Date.fromYMD(y,m,1))

});
}).flatten());

print(monthlyImages)

How do I do a decadal composite (every 10 days)?


texts: Finding latitude and longitude limits for a sub-region given its name

texts:
In order to extract climate data from Copernicus, there's the Sub-region extraction section that allows one to extract data from a region within specific coordinate limits.
I am wondering how should one do to find the latitude and longitude limits for a specific sub-region, given the name of that region.
Taking the Iberian Peninsula as an example, what should one do in order to find its North, South, East and West limits?


texts: How can I use UpdateCursor in AGOL layers?

texts:
I'm all new to Python, and 'Im stuck with a problem:

Is it possible to use UpdateCursor on Feature Layers in AGOL?

I am trying something like this:

from arcgis.gis import GIS

gis = GIS("https://arcgis.com", "xxxxxx")

import arcpy

fs=gis.content.search(query="title:xxx, owner:xxxx", max_items=5000,)

kort = fs.layer[0]

fields = ["Fabrikat", "Entreprenør"]

with arcpy.da.UpdateCursor(kort, fields) as cursor:

    for row in cursor:
        if row[0] == "ABB":
               row[1] = "3"
        cursor.updateRow(row)

but its popping up with this error:

  RuntimeError                              Traceback (most recent call
  last)  in
        2 lag = query_result1
        3 fields = ["Fabrikat", "Entreprenør"]
  ----> 4 with arcpy.da.UpdateCursor(lag, fields) as cursor:
        5     for row in cursor:
        6         if row[0] == "ABB":

RuntimeError: 'in_table' is not a table or a featureclass

How can I make the cursor "eat" the layer?


texts: Why ExtractMultiValuesToPoints not working for Point file?

texts:
I'm struggling to extract depths from a bathymetry raster to a Point file with multiple features (not Multipoint, as I previously mis-stated). In trying to use ExtractMultiValuesToPoints(), I noticed it worked on feature classes with one point, but not those with multiple points. Here's an example I've run on two test files.

import arcpy
import arcpy.sa as SA

bathy = [filepath]
testgrid1 = [filepath] (one point)
testgrid4 = [filepath] (four points)

SA.ExtractMultiValuesToPoints(testgrid1, [[bathy, "depth"]], "BILINEAR")
SA.ExtractMultiValuesToPoints(testgrid4, [[bathy, "depth"]], "BILINEAR")

Results: "depth" columns are added to both files. The correct depth inserted into the depth column of 'testgrid1'; but  values are inserted into the depth column of 'testgrid4'.

I'm certain that there's data in the bathymetry file at the location of the points, so I can't fathom why this isn't working.


texts: Combined Hstore key/value and spatial query too slow to handle bigger OSM extracts

texts:
I'm trying to compute some statistics for OSM data using PostgreSQL 9.3.5 and PostGIS 2.1.4. I started with a small bavaria extract which I downloaded from Geofabrik. The db schema is the normal API 0.6 schema, the data was imported via the dump approach into Postgres (using the pgsnapshot_schema_0.6*.sql scripts which come with osmosis). ANALYZE VACUUM was also performed.

The only custom made thing I'm using is a polygon table which contains multipolygons for all administrative boundary relations. The polygon geometry was not simplified in any way.

What I'm now trying to achieve is counting all nodes which are inside of the admin=6 boundaries of bavaria. Here is my SQL query:

SELECT relpoly.id, count(node)
FROM bavaria.relpolygons relpoly, bavaria.nodes node
WHERE relpoly.tags @> '"boundary"=>"administrative","admin_level"=>"6"'::hstore
AND ST_Intersects(relpoly.geom, node.geom)
GROUP BY relpoly.id;

The runtime of this query is terrible because Postgres is doing a nested loop join and scans over all nodes for every admin=6 boundary. FYI, bavaria is divided into 98 admin=6 polygons and there are about 30 million nodes in the bavaria extract.

Is it possible to avoid this sub-optimal query execution and to tell Postgres that it should scan all nodes only once (e.g., by incrementing a counter for the corresponding polygon in the result set or by using hints)?

Edit:

1) a spatial index exists on the bavaria nodes:

CREATE INDEX idx_nodes_geom ON bavaria.nodes USING gist (geom);

2) the query plan look like this:

HashAggregate  (cost=284908.49..284908.75 rows=26 width=103)
  ->  Nested Loop  (cost=111.27..283900.80 rows=201537 width=103)
        ->  Bitmap Heap Scan on relpolygons relpoly  (cost=4.48..102.29 rows=26 width=5886)
              Recheck Cond: (tags @> '"boundary"=>"administrative", "admin_level"=>"6"'::hstore)
              ->  Bitmap Index Scan on relpolygons_geom_tags  (cost=0.00..4.47 rows=26 width=0)
                    Index Cond: (tags @> '"boundary"=>"administrative", "admin_level"=>"6"'::hstore)
        ->  Bitmap Heap Scan on nodes node  (cost=106.79..10905.50 rows=983 width=127)
              Recheck Cond: (relpoly.geom && geom)
              Filter: _st_intersects(relpoly.geom, geom)
              ->  Bitmap Index Scan on idx_nodes_geom  (cost=0.00..106.55 rows=2950 width=0)
                    Index Cond: (relpoly.geom && geom)

3)

I created the following two indexes, but the query plan (and runtime) did not change

CREATE INDEX relpolygons_tags_boundary on bavaria.relpolygons( (tags->'boundary') );
CREATE INDEX relpolygons_tags_admin on bavaria.relpolygons( (tags->'admin_level') );
ANALYZE bavaria.relpolygons;


texts: Python Script Performance as an ArcGIS Tool Versus Stand-Alone

texts:
Has anyone studied the difference in running a Python script in ArcToolbox versus as a stand-alone script? I had to write a quick-and-dirty script to convert a set of RGB images to single band by extracting band 1. As a stand-alone script reading and writing to my PC it processes 1000 identically-sized images in about 350 seconds. Running the same script from ArcToolbox takes about 1250 seconds.

import arcpy
import csv
from os import path

arcpy.env.workspace = in_folder
image_list = arcpy.ListRasters()

#Create a CSV file for timing output
    with open(outfile, 'wb') as c:
        cw = csv.writer(c)
        cw.writerow(['tile_name', 'finish_time'])

        #Start the timer at 0
        start_time = time.clock()

        for image in image_list:
            #Extract band 1 to create a new single-band raster
            arcpy.CopyRaster_management(path.join(image, 'Band_1'), path.join(out_folder, image))
            cw.writerow([image, time.clock()])

I added some code to track when each tile finishes processing, and export the results as a CSV. Converting the finish time to processing time occurs in Excel. Graphing the results, the processing time is roughly the same for each tile as a script, but processing time increases linearly when run as an ArcGIS Tool.

If the data reads and writes are to a network device, the increase appears to be exponential.

I'm not looking for alternate ways to accomplish this particular task. I want to understand why the performance of this script degrades over time when run as an ArcGIS tool, but not as a stand-alone script. I have noticed this behavior with other scripts as well.


texts: Importing paths for QGIS 3 standalone scripts

texts:
For QGIS 2.x, I used the batch file provided in the answer for this post in the OSGeo4W Shell to run standalone scripts: Problem with import qgis.core when writing a stand-alone PyQGIS script

The batch file for QGIS 2.18.17-ltr looks like this:

set PYTHONPATH=C:\Program Files\QGIS 2.18\apps\qgis-ltr\python
set PATH=C:\Program Files\QGIS 2.18\apps\qgis-ltr\bin;%PATH%

Now after installing QGIS 3.0, I would also like to run some standalone scripts. So I tweaked the above batch file:

set PYTHONPATH=C:\Program Files\QGIS 3.0\apps\qgis\python\
set PATH=C:\Program Files\QGIS 3.0\apps\qgis\bin;%PATH%

However, when I load this batch file into the shell and use from qgis.core import *, I receive the following error message:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Program Files\QGIS 3.0\apps\qgis\python\qgis\__init__.py", line 26, in <module>
    from builtins import zip
ImportError: No module named builtins

How can I set the paths correctly?


texts:
I want to check a vector file if there is a column "Source".

If yes, nothing to do. If there is no column "Source" it should be created and filled up with '0'.

Is this possible to do via the modeler?


texts: OpenLayers sending flipped BBOX coordinates when using custom projection

texts:
Following the documentation found here, I've been trying to create an OpenLayers map view using a custom projection, EPSG:4674 to be more specific. My current code looks like this:
import proj4 from 'proj4';
import {register} from 'ol/proj/proj4';
import {get as getProjection} from 'ol/proj';

proj4.defs('EPSG:4674','GEOGCS["SIRGAS 2000",DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6674"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4674"]]');
register(proj4);
const SIRGAS2000 = getProjection('EPSG:4674');

const map = new Map({
  view: new View({
    projection: SIRGAS2000
...

and those are my layers, in my GeoServer instance:

The problem happens whenever I try to load any WMS source using the following code:
import { Tile as TileLayer } from "ol/layer";
import { XYZ as XYZSource, TileWMS as TileWMSSource } from "ol/source";

const estadoLayer = new TileLayer({
  title: "Estado",
  source: new TileWMSSource({
    url: GEOSERVER_WMS_URL,
    params: { LAYERS: "postgis:estado", TILED: true },
    projection: SIRGAS2000,
  }),
});

OpenLayers sends reversed BBOX coordinates to my GeoServer, as you can see here:

The funny thing is: If I remove the projection attribute from TileWMS constructor and set View's projection to "EPSG:4326", everything works fine:

What could be causing this issue? Considering that EPSG:4674 should, in theory, be compatible with EPSG:4326.


texts: QGIS Layer Tree API crashing when removing layer

texts:
I'd like to change the order of layers in my layers panel using Python in the Python Console (this is a for a script not a plugin).

Following a blog post here on the Layer Tree API for QGIS, the following code causea QGIS to crash:

root = QgsProject.instance().layerTreeRoot()
testClone = root.children()[0]
root.insertChildNode(1, testClone)
myLayer = root.children()[0]
root.removeChildNode(myLayer)

I've tested this on QGIS 3.4.15 and 3.8.3. Both crash. Manually removing the original layer also causes a crash.


texts: Remove all but specific highways in osmosis

texts:
I want to display entire country OSM data in my application.
The problem my application don't know how to deal large maps.
Instead of chunking the map into 5-10 small chunks, I wish to display only the following items:

highway=trunk, highway=trunk_link, highway=primary,highway=primary_link, highway=secondary, highway=secondary_link,highway=tertiary,highway=tertiary_link ,highway=motorway, highway=motorway_link

All roads but residentials. I don't want to show rails, polygon shapes (area, landuse, parks.. etc..) Only roads!

But after running this command:

osmosis.bat --read-xml "C:\belgium-latest.osm" --tag-filter accept-nodes highway=trunk, highway=trunk_link, highway=primary,highway=primary_link, highway=secondary, highway=secondary_link,highway=tertiary,highway=tertiary_link ,highway=motorway, highway=motorway_link --write-xml "C:\belgium-latest-highwayOnly.osm"

I still have the other shapes (not only highways) and the osm file is too big.. it shrunk only 25%~.. i hoped it to be smaller.

I am doing something wrong ? should i add objects to reject-nodes attribute ?


texts:
We are using the Collector for ArcGIS app for numerous inspections. I've been tasked to manage the GIS portion of it.

Some locations lack data service so the inspectors are using the offline mode. I have been asked to add a feature service to the map after the fact and cannot seem to get the app to pull the new feature service into the offline mode.

Do the users have to re-download the map to access this new editable feature service?

I was under the impression that syncing it should pull new features as well.


texts: export WMTS tiles as MBTiles with reprojection

texts:
I need to export tiles from a WMTS server to the MBTiles format. The server is not public, but I have access to it and permission to export a specific portion of the tiles. The problem I have seems to come from the different projections:

The WMTS server uses UTM zone 38N (EPSG:23038)
The exported MBTiles needs to be in WGS84 spherical mercator (EPSG:3857)

What's the easiest approach to achieve this, and with which tools?

I already tried loading the WMTS source in QGIS and using QTiles. This works, but the resulting tiles are rendered poorly, I guess due to a sub-optimal reprojection algorithm.


texts: Uniquely identify feature classes & compare them for equality using ArcObjects?

texts:
I want to implement EqualityComparer<T> for feature classes so that I can stuff them into a Dictionary or HashSet.

For motivation, say that I am looping over many IFeature objects which belong to different feature classes, and I want to set field values for each. FindField is inefficient to perform inside of a loop. If I have a Dictionary<IFeatureClass, Dictionary<string, int>> to give me the index of a given field, this operation would be much faster.

So, how does one uniquely identify a feature class object? The AliasName property is not unique. The IName object should be, but I'm not sure how to compare them, or how to hash them. Is IDataset.Name unique to each feature class? If it is, it must include the full data path, right?


texts: In DIVA-GIS, how do I set a .grd to show as containing continuous or categorical data?

texts:
I am running both DOMAIN and Maxent models on a set of raster data. In Maxent, I am able to set whether each dataset contains either categorical or continuous data. Is there a way to do this in DIVA-GIS before I run DOMAIN?

I have found this question in a few places, without any answers. I have read one thesis where the variable was set by using a hex editor, but that paper is over 10 years old.

Bottom line, I want to be sure the model is set up properly so I get the best results.


texts: Is there any existing postgis package for Ubuntu 16.04 (xenial) compiled against GDAL 2.0+?

texts:
I stumbled upon this error while using QGIS (2.18.13) to display a PostGIS raster and read in comment #14 on the same issue, that it is working when using a newer GDAL version.

I have GDAL 1.11.3 used by the DB, which I installed using this instruction (PostgreSQL 9.6.5 & PostGIS 2.3.3). I am curious, if anyone has a PostGIS package for Ubuntu 16.04 compiled against a newer GDAL version. Or do I need to compile it from scratch?


texts: Spacing points automatically instead of by hand using ArcGIS Desktop?

texts:
Is there a way to put points out at 10m apart instead of doing it all by hand?

My boss thinks it can be done in Excel but the only way I know how is to enter in each lat, long then import it that way.


texts: Geoprocessing Tools Create Empty Shapefiles in QGIS?

texts:
I seem to have great difficulty using the geoprocessing tools in QGIS. Sometimes they work fine, sometime no matter what I do I get a null result, even though I can see I should be getting a selection. I am mostly using overlap and clip, and I have made sure the two (polygon) shapefiles use the same CRS and that if a selection is used that it overlaps.

I must make it clear that this does seem to be an intermittent problem, but restarts don't cure it.

Is there anything else I should be looking for?


texts: Unselecting Individual Selected Features using ArcObjects?

texts:
In my Add-In built using .NET ArcObjects/C# with Arc 10.0 I would like to be able to programmatically unselect a subset of all selected features gathered through a reference to  FocusMap.FeatureSelection. I can pretty easily holistically unselect all of the features using the Clear() method, but I need to do this on a more targeted basis. For example, I display a list of all selected features (obtained by casting the FeatureSelection ISelection to an IEnumFeature) to the user and ideally would like to be able to provide some management (read: removal) of the selected features to the user.

So, in a nutshell: Is there a way to remove just a certain feature or features from an entire selection in ArcMap using ArcObjects and if so, what is the preferred method?


texts: Adding information from one polygon to the attribute table of another

texts:
I would like for the watershed (blue) attribute table of this map to have a column that indicates what protected area (red) it overlaps with. (As you can see, multiple watersheds touch each protected area).
I have searched around and am unable to figure out how to do this using ArcMap 10.1.  Any suggestions or links to solutions would be awesome!


texts: How to create multiple points inside polygon in QGIS

texts: How do I generate five random points within each of three polygons with a minimum distance of 100m? I am using QGIS 2.18.4 and 'Random Points Inside Polygons (Fixed)' Sampling Strategy 'Point Count' 5 and min dist 100. I am currently only returning a single random point in each polygon rather than 5. There is plenty of space in each polygon to locate 5 points with >100m between. Additional information: Project and Polygon are EPSG:28356, unit set to meters.


texts: Creating landslide susceptibility map using logistic regression and ArcGIS Desktop?

texts:
My project is to build a landslide susceptibility map using logistics regression. I have generated the following maps so far:

slope map, hill shade, curvature, elevation, road buffers, river buffers of 100 - 500m , Aspect, land use land cover, ndvi.

All of these maps after production do not have attribute tables. please, I need a detailed explanation on how to use some of this data produced to produce a susceptibility map using logistic regression.

I'm using ArcGIS Desktop 10.0.


texts: Arcpy: Names and distances of nearest intersecting roads in layer 1 to end points of a road segment from layer 2

texts:
As shown in the image below, I'd like to return a new field in layer 2 that for each segment stores a string:"Road C, 1 km from Road A to 0.5km From Road B".

Calculation of the distances could be considered optional; however, it must be able to return the nearest road name to each of the segments vertices. Names are just an attribute of Road C.

If there is no nearest road on the end side (Dead end) then the script should return None.

.

What I have been able to do so far is generate a nearest table that gives me a nearest FeatureID. I join that back to the Layer 1 and return a name; however, that doesn't give me a distance or the nearest road on the other end.

Note that all of these roads are part of the linear reference system

EDIT 1:

the solution cannot pick an intersecting road that's within the "green" road's boundaries.


texts: ArcGIS Server Sporadically Exports Blank Image / renders blank image, for a published map service

texts:
For several months our map services in our ArcGIS Server 10.5.1 setup occasionally fail to draw/render when browsing to their export REST endpoints.

Our setup is in the cloud and we use cluster setup. We also are using SDE. We have tested and verified that the issue is not because of any load or any service configuration (such as recycling time, etc.) or network issue or firewalls, etc. It just happens sporadically. We haven't been able to pinpoint the cause of this behaviour. All thise using HTTPS.

Issue:

When browsing the REST endpoint of a published map service, a blank image is rendered sporadically. In other words, image gets displayed correctly sometimes and sometimes image is blank. So on the client side web interface, when user browse the map, a map service will render correctly and as soon as user pans, the map service will not display. Again when user pans the map, map service will render correctly. This happens at anytime (i.e. nothing fixed or no pattern observed).

We have observed that this happens to particular services and haven't seen such wierd behaviour for other map services.

Investigation carried so far:

Restarting the bad map service in Server Manager or ArcCatalog.
Restarting the ArcGIS Server Windows service.
Restrting the cloud servers VM itself.
Trying various service configurations such as recycling intervals, etc.
Load testing and automatic querying the map service with many requests.
Creating MXD from scratch and publishing it.
Verifying MXD by using simple symbology or not using complex definition queries or complex labelling, etc.
Checking data integrity.
Verifying for reliable network.
Verifying if one map service is failing to rendering then other services are working as expected.
Increasing heap sizes.
Checking AGS logs (no errors encountered).

Has anyone heard of or seen something similar to this? And has anyone resolved it or found the cause of this?


texts: Export GeoTIFF in lat long not map units?

texts:
I'm trying to export a map as raster (GeoTIFF with world file) from print composer. The world file is always coming with Map units because I'm using Google Satellite from OpenLayers web plugin. I tried to change the CRS from Projected coordinate (EPSG:900913) to Geographic coordinate (edited : WGS84 EPSG:4326) in order to get lat/lng in the world file and the world file come as expected. But after changing to EPSG:4326 the position of some features were not correct. So I assume the transformation I'm doing is wrong.

The objective is to use Google basemap, add my layers and labels, then save everything as GeoTIFF to be used in other software. but that software is accepting only GeoTIFF with lat/long info in the world file not projections as map units.

Is there other way to do this conversion?


texts: Getting large amounts of bounding boxes around cities and natural land coverages

texts:
I am an astronomy student working on a research project in which I want to get large amounts of bounding boxes around cities and natural land coverages. In QGIS, I am using a dataset of Europe that gives the different land coverages (types of forest, cities, beaches, you name it) and of course the option to turn those on or off in display.

I would very much like to automate the process of selecting a certain land coverage type, and retrieving a whole bunch of BBox values (maybe within a certain subset of the entire dataset) to then use in my research. I would very much like to ask if any of you good give me any headers about what route to take to achieve this.

I am used to working with Python, and I think that writing a Python script is maybe a way to go about this? I am completely new to working with GIS, and hope this is the right place to ask.


texts: Is there an easy way to search for a road name in a polyline layer in QGIS?  I'm looking for something like what ArcMap 10.1 has, just a basic search box where I can type in a name and have it give me a list of matches, then have the option to zoom to it, or highlight it, or something


texts:
Suppose I have a field titled NAME1.

Examples of what is in the field appear like this:

Fox Lake 1
WABASCA 166
TSU TUE 196G
YELLOW QUILL I.R. 90-8
POUNDMAKER INDIAN RESERVE NO. 114-17
CROSS LAKE 19X02
BEARDY'S & OKEMASIS I.R. 96 & 97A

How using python or vb script would I extract a string only after the first digit into a new field titled Name_Num?

Basically, the Name_Num field would have to appear like this:

1
166
196G
90-8
114-17
19X02
96 & 97A

Mind you, I need the letter attached to the numbers. I can't just extract digits only. But anything before the first number must go. The numbers vary in character length but always appear at the end of my string.

I rarely use python and am not familiar with when something goes into the pre-code block or doesn't.

These are the results after inputting the first reply which looked great in my opinion but didn't work unfortunately.


texts: Pyqgis - Removing a layer from the composer legend

texts:
I am developing a script to print different map views of a QGIS project.
I want to print all the vector layers in the canvas but I don´t want one of them to appear in the legend. I can`t figure out the way to do it.


texts: loading multiple data in a delimited text file

texts:
I am trying to map out a pipeline, identifying where the pipeline is covered by more than 0.6m of sediment and where the pipeline is not covered by 0.6m of sediment. I have gathered this information in excel and saved it as a csv. file. (seen below) I also have the same excel sheet expressing DOC (depth of coverage) being greater than 0.6m.

Note, I have got over 1000 data points all with a specified easting and northing. Using QGIS I add layer --> Add Delimited Text Layer and upload the document (Shown below).

When I add the layer, It doesn't show on my map even when I right click the layer and zoom to layer.

What I am trying to achieve is to plot all 1000 data points according to their Easting's and Northing (in a line), where I can change them from point data to line data then I can add rules to it so I can identify where the areas lower than 0.6m are and where they are greater than 0.6m. Does anyone know a way of doing this?


texts:
I'm creating a style in Mapbox which is supposed to contain hundreds of shapefiles for a total size (shp + dbf) of several GBs. Due to the limit of 15 sources per style, I've combined the smaller shapefiles into shapefiles of about 250 MB, and styled the data from the different shapefiles accordingly to the name of their original shapefiles.
Unfortunately, this is not enough for representing all my data within the style because I exceed the limit of 15 sources x 260 MB each.

Am I doing the above correctly or am I missing something that would allow myself to use more sources?

EDIT: I'm experimenting with Tippecanoe and MBtiles. Hopefully it will not oversimplify my shapes.


texts:
Firstly, sorry for having to post this question, I feel like I've been going in circles trying to make sense of similar questions. I'm quite new to QGIS so please do re-direct this if needed.

Why I want to do I want to crop the NLCD data raster, e.g. Raster -> Extraction -> Clip Raster by Mask Layer so I can do some very basic counting of land-types on a per-county basis. This in itself seems to cause issues as the Raster data is just so big. So if there is some other method that I am unaware, please let me know. Otherwise, please read on.

Background: Like many that have come before me, I have been using the NLCD Land Cover CONUS Data set for the US. This provides a pretty large set of files (~20 GB per year), of which a erdas .img file is easily created in QGIS (drag and drop). In the metadata file, this has the projection:

PROJCS["Albers_Conical_Equal_Area",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-96.0],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_origin",23.0],UNIT["Meter",1.0]]

Next to this, we have the US 2017 Census County Shapefiles. Again this is easily entered into QGIS with a drag and drop. This has the projection:

GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

Problem Summary: Default QGIS CRS is 4326. I load in the NLCD data, set the NLCD layer CRS to EPSG:102003 (Albers). I then set the project CRS to the same. Then loading in the Shapefiles for the US County data, I get the image below. Showing that the layers are offset.

I feel like the issue could relate to the problems as described here and here but I'm not exactly sure what is going wrong in my process. Further, my data is so big that any warping seems to take forever on my laptop.

Am I just missing something obvious?


texts:
I am trying to label Census Block Groups and Tract ID for the block group feature layer. I would like to label them using different color halos. Is there an arcade expression that will allow me to change the color of the halos for each label in a feature?
So far, I have been able to change the color of the text for the labels but not the color of the halo.
"<CLR red = '0' green = '197' blue = '255'>" +
'Tract ' +
"<BSE>" + $feature.TRACT + "</BSE>" +
"</CLR>" +

TextFormatting.NewLine +

"<CLR red = '223' green = '115' blue = '255'>" +
"BG - " + $feature.BLKGRP +
"</CLR>"

From the image below, I basically want to invert the colors of the halo and the text.


texts: varying size of symbol depending on number of datasets

texts: i am trying to map archaeological finds with qgis. some of them have the same coordinates since they are from the same site or very close to each other. is it possible to automatically display the number of finds in a close range through a varying size of the symbol? so, the more datapairs i have in the same spot, the bigger the symbol gets? otherwise 10 finds just appear as one in the map because they are kind of on top of each other...


texts: Update Data Driven Pages Index layer via Python

texts:
I am trying to find a way to update the index layer for data driven pages through arcpy. I understand that to access DDP from arcpy.mapping you need to be referencing a MXD that has DDP enabled. My problem is that I am trying to create DDP for a feature class I just created from the XML to Geodatabase import tool and use that as my index layer instead of the scratch layer I created to enable DDP. The MXD contains a layout I would ultimately like to use when exporting my DDP. Here is my code so far.

# Import arcpy module
import arcpy, os
from arcpy import env

env.overwriteOutput = True
env.workspace = r"C:\Users\spadgett\Desktop\ArcScrap\Tyler_TEST.gdb"

# Script parameter
xml = arcpy.GetParameterAsText(0)

# Hardcoded GDB for storing xml
Tyler_TEST_gdb = "C:\\Users\\spadgett\\Desktop\\ArcScrap\\Tyler_TEST.gdb"

# execute xml to geodatabase tool
arcpy.AddMessage("Executing XML to Geodatabase")
arcpy.ImportXMLWorkspaceDocument_management(Tyler_TEST_gdb, xml, "DATA", "")

#create feature layer for points
FCpoints = "Control_Point_Name"
arcpy.AddField_management(FCpoints, "SCALE", "SHORT")
arcpy.CalculateField_management(FCpoints, "SCALE", "1200", "PYTHON")
arcpy.AddMessage("Points feature class added with SCALE field")

arcpy.MakeFeatureLayer_management(FCpoints,"Points2")

# Set environment and MXD
mxd = arcpy.mapping.MapDocument(r"C:\Users\spadgett\Desktop\ArcScrap\Tyler Sample\Tyler_Test.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]

layer = arcpy.mapping.Layer("Points2")
arcpy.mapping.AddLayer(df, layer, "TOP")

RefPoints = r"C:\Users\spadgett\Desktop\ArcScrap\LYR\Points.lyr"
updateLayer = arcpy.mapping.ListLayers(mxd)[0]
arcpy.ApplySymbologyFromLayer_management (updateLayer, RefPoints)
mxd.save()

updateLayerOld = arcpy.mapping.ListLayers(mxd)[1]

arcpy.mapping.UpdateLayer(df, updateLayerOld, updateLayer, False)
mxd.dataDrivenPages.refresh()
mxd.save()

#Create JPG of DDP for Points name
##mxd = arcpy.mapping.MapDocument("CURRENT")
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
    mxd.dataDrivenPages.currentPageID = pageNum
    pageName = mxd.dataDrivenPages.pageRow.Name
    arcpy.mapping.ExportToJPEG(mxd, r"C:\Users\spadgett\Desktop\ArcScrap\Tyler Sample\Images\AGM_" + str(pageName) + ".jpg")

Running this code ends up with an error at the UpdateLayer line. Switching the boolean to True just changes the symbology, just as the help document says.

Any ideas?

Error is

Traceback (most recent call last):
  File "C:\Users\spadgett\Desktop\ArcScrap\Tyler Sample\python\AGM_Report.py", line 59, in <module>
    arcpy.mapping.UpdateLayer(df, updateLayerOld, updateLayer, False)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\mapping.py", line 1876, in UpdateLayer
    tc.ReplaceLayer(tl, rl)
ValueError: DataFrameObject: Unexpected error


texts: Multiple GeoServers with Integrated GWC for tile seeding

texts:
I am trying to run a pretty large seeding task, so have setup multiple tomcat 8 services running on different ports and deployed the latest geoserver WAR using a single data directory using the environment variable.

However when I then use the integrated GeoWebCache and seed an area it seems they overwrite each other even though the seeding type to SEED (rather than reseed).

Having read some more the of the docs there seems to be a Tile Locking Mechanism called NIO that can used.

NIO File locking: Uses Java New IO File Locks suitable for use in a clustered environment (with multiple GeoServers sharing the same data directory).

So on all my GeoServers I changed the configuration and re run a sample area and it still seems the different GeoServers are running independently.

In my GWC folder I do get a new lockfiles folder but nothing appears in it.

Maybe it is not working as I have set a different folder for the GWC tiles in the web.xml file rather than the shared data directory I used for the styles, workspaces and stores mentioned above.

Can anyone confirm or provide a way to get this to work?


texts:
I use QGIS 2.12 64-Bit Version on Win7.
I have a problem with red matrix wich show up around every label. I want to delete them or make them invisible. It looks like in the following screenshot

Does somebody know how to solve this?


texts: Split a polygon into smaller parts without using Split

texts:
Is there an alternative to split in ArcMap 10.4, if you don't have the advanced license?

What I want to do is split a polygon according to another polygon. I want to preserve the extent of the original polygon but have it be divided into smaller parts.

In the image below, the underlying polygon (looks like a raster but is a polygon) should be divided according to the borders (white lines) of the overlying polygon.

With the Clip tool, you extract parts of the polygon instead of just splitting it, so that doesn't seem like an option.


texts: what is the main idea of new (since osm2po ver 4.7.7) MultiLevelGridRouter

texts: in the new osm2po v4.7.7 release notes I see that there is new MlgRouter (MultiLevelGridRouter) "which is an at least five times faster DefaultRouter". Could someone clarify why it's so fast (what assumptions/concept cause such dramatic speedup? what is the main idea of Multi Level Grid?) and for what cases it's applicable/recommended to use?


texts:
I'm trying to find information about the possibilities to map an arc (or a circle to be specific) as a straight line. Imagine a circle mapped onto a country or continent, that represents a path with a starting point 'A'. The circle has a small buffer. I'm trying to map this path into a straight line starting and ending with point 'A', so that I end up with, instead of a square map containing the circle, a very long (length of the path) and thin (height of the buffer) rectangle shaped map.

I figure that I will need to create a custom projection of some sort, and possibly splitting the circle into two or more arcs. But I have a hard time finding a grip on the subject, and think I might be overlooking a simple solution.


texts:
I'm using ArcScan in ArcMap version 10.8 with a raster I exported from GEBCO.  I have the data classified as black for anything less than zero (ocean/bathymetry) and white for anything greater than zero (land).
When I open the Raster Painting toolbar (via the Raster Cleanup dropdown in ArcScan) and select the raster Brush, it allows me to "paint" and change from foreground to background color, but when I release my mouse click (ie "lifting the paintbrush") the changes never persist.
How do I get the Raster Brush "paint" to persist and save?


texts: TableToTable_conversion using PostgreSQL and ArcPy

texts:
Currently, I am trying to convert a geodb to a table in my PostgreSQL db with arcpy as such. I've tried both with the parcels table being created and the parcels table not being created

arcpy.TableToTable_conversion("C:/Temp/geodb.gdb\parcels", "C:/../postgresconn.sde", "parcels")

I know my ../postgresconn.sde works because I can do other things such as FeatureClassToFeatureClass_conversion via...

arcpy.FeatureClassToFeatureClass_conversion("C:/../postgresconn.sde\database.public.area", "C:/Temp/geodb.gdb", "area", "", "", "")

The problem is in the TableToTable method, is that I cannot use database.public.area (dots are not allowed in that argument string)

I have set my searchpath on the PostgreSQL sql side to point to database.public hoping that it would just reference that schema.

The error I am getting ...

arcgisscripting.ExecuteError: ERROR 000210: Cannot create output C:\...\postgresconn.sde\area

How can I make the method call TableToTable work in this scenario?


texts: Chemical Spill on Land (ArcGIS, ArcMap 10.4.1, Hydrology Tools)

texts:
I'm a very new ArcGIS user currently working on my first project regarding chemical spills on land.

The Background:

Software that I have to use is ArcMap.

Currently, I am trying to use the pour-point tool to delineate the flow path from a proposed spill/leak point of a chemical bearing pipeline.

So far I have used only the Hydrology Toolset and this is what I did in sequence...

Filling the LiDAR Data
Generate Flow Direction
Generate Flow Accumulation from the Flow Direction
Delineate the drainage Basin and also turned it to polygon features
to outline the flow accumulation passages.
I created a new point shapefile and assume it as the "proposed
leaking point" (Please view the hand drawn diagram below)

The problem:

Firstly I measure by estimating the distance of the "proposed leaking point" to the nearest flow accumulation line (12m). Then I run the snap-pour point tool with a little addition to the measurement (snap-distance: 15m). After it transforms the point shapefile of the "proposed leaking point" as a pour-point (as shown in blue in the diagram below) and snapped it to the nearest highest flow accumulation point... I then use the watershed tool..

However... This only delineates the whole stream (watershed) which contributes to the pour-point. This is not what I am looking for as it is not taking into account the "proposed leaking point".

The Question:

How do I get the delineation or the flow line of the chemical from the "proposed leaking point"?

ArcMap is not working so I had to draw this.

Edit: Additional Images

 Image above shows the pour-point that creates a watershed of the contributing streams.

 Image above shows where I would like the flow line to be if there is a leak at that point of the pipeline.


texts: Where command in raster calculator Arcgis

texts:
I have two rasters that I wish to combine to give a third raster that describes spatially where multiple conditions of both rasters are met. For example, say my two input rasters are raster A and raster B. Values like this:

RASTER A:

   X1 X2 X3
Y1 A1 A2 A2
Y2 A1 A1 A1
Y3 A1 A1 A2

RASTER B:

   X1 X2 X3
Y1 B1 B1 B1
Y2 B1 B1 B1
Y3 B1 B1 B1

I want to create a third raster with different values that occur where conditions from both are met, see the below "psuedo-code" describing what I want:

WHERE (RASTERA = A1 AND RASTERB = B1) THEN ALPHA
WHERE (RASTERA = A2 AND RASTERB = B1) THEN BETA

Giving my final raster as below:

   X1    X2    X3
Y1 ALPHA BETA  BETA
Y2 ALPHA ALPHA ALPHA
Y3 ALPHA ALPHA BETA

How can I accomplish this, preferably using raster calculator?

P.S. In my real data, I have numbers that I need to be in certain ranges. As such I could not simply concatenate the raster values (e.g. A1, and B1 to give A1B1) and give an alias fitting to that.


texts:
I have a React JS component that adds a Layer to the Basemap:

 componentDidMount() {
    const { client } = this.props;

    // Craete source, styles and layer with the given props
    const cartoSource = new carto.source.Dataset("dataset");
    const cartoCSS = new carto.style.CartoCSS(style);
    const layer = new carto.layer.Layer(cartoSource, cartoCSS, {
      featureClickColumns: [
        "geoid2"
      ]
    });

    // Add them to the client and to the map
    client.addLayer(layer);
    client.getLeafletLayer().addTo(this.context.map);

I would like to display the geoID when I hover over the polygon, and also highlight it. Having a hard time understanding how this could be done in React.


texts: Using Union in ArcPy for ArcGIS Pro results in two fields with same name?

texts:
I just performed a union in my Python script between 2 features that both had a field attribute named "FLD_ZONE". My hope was that one of those attributes (for the second feature, I assumed) would be renamed "FLD_ZONE_1" in the union, but this was not the case.

Is there a way using Python to rename that second "FLD_ZONE" field name to "FLD_ZONE_1" or something different?

I know I can do that manually, but I'm not sure how that would work in a script.


texts: Low-res coastline of world landmass, in xml format?

texts:
Can anyone direct me towards an xml file, in typical OSM format, for the low-res coastal Ways of every land-mass on Earth?

I need something easily renderable at low zoom levels, avoiding the superfluous detail of full-res data.


texts: Output extend for raster calculator (setnull set)

texts:
We are working on a project to find least cost path for marine shipping based on some criteria. In a phase when iterating raster datasets and using raster calculation to use setnull function to exclude some criteria in a final raster. Problem with the output is that the extend is based on a smaller raster (one dataset is for the whole world and the second one just for arctic regions which the output is clipped). We want the output to appear on a worldwide scale.
I have tried to change everything in geoprocessing environments but still getting output just for the area from a smaller raster.


texts: Authentication in GeoServer with OpenLayers2

texts:
Is there a way to authenticate in GeoServer with a external webpage (made with asp.net and OpenLayers)? I have restricted the access data only by certain users, and I don't want to authenticate this way:

username:password@geoserverurl

for security reasons (the username and password are visible in the webpage code)


texts: Changing of unknown CRS of the shapefile in QGIS

texts:
I'm trying to convert "unknown" CRS of the shp to the WGS84 UTM zone N34. I've tried to do it with saving of shp in different CRS (Save as) but I did not succeeded. I've checked the situation with open layers and using of on-the-fly transformation and matching is far from perfect (500-600 m difference).

Unknown CRS has been showed in QGIS as USER 100001 with the following parameters:
 * Generated CRS (+proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +units=m +no_defs)

It should be a Gaus Kruger projection (Bessel 1841) for Serbia. I've tried also with changing of shp and open layer to MGI 1901/Balkans zone 7 and MGI/3-degree Gauss zone 7 but the situation is the same.

Note: Shapefiles should be very precise as they were obtained by the Public Enterprise. Unfortunately, they're not interest in delivery of the shapefile in a different CRS.

I would be very greatful if someone could help me regarding this issue.


texts: Code to load TIFF image on openlayer.layer.wms using MapServer MS4W

texts:
Where did I get mess up?

Here is my code inside the mapfile.map

LAYER
    NAME "layername"
    STATUS OFF
    TYPE RASTER
    DATA "data/sample.tif"
END

The code inside explorer.js from where I declare the openlayers.layer.wms

    var showlayer = new OpenLayers.Layer.WMS(
        "my Layer Title",
         my_url,
        {
            layers: 'layername',
            'format':'image/tiff'

        },{
            singleTile: false,
            maxExtent: my_MapExtent,
            isBaseLayer: false

   })

   map.addLayer(showlayer);

When I open it in the browser here the error show's:

msWMSLoadGetMapParams(): WMS server error. Invalid layer(s)

any help or suggestion for this one?


texts: Converting point cloud (.txt) to shapefile with ArcPy?

texts:
I have lots of .txt  files (point cloud) which are raw data acquired from from LIDAR and containing XYZI values. I want to convert them to shapefile.

Is there a way I can do that instead of adding XY data first and exporting to shapefile one after the other?