texts:
I modified my repository like this:
deb https://qgis.org/ubuntugis xenial main
deb-src https://qgis.org/ubuntugis xenial main

Then:
sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable

Finally:
wget -O - https://qgis.org/downloads/qgis-2017.gpg.key | gpg --import
gpg --fingerprint CAEB3DC3BDF7FB45^C

I got a lot of broken packages:
Obj:1 http://mirror.cedia.org.ec/trisquel.packages flidas InRelease
Obj:2 http://mirror.cedia.org.ec/trisquel.packages flidas-security InRelease
Obj:3 http://mirror.cedia.org.ec/trisquel.packages flidas-updates InRelease
Obj:4 https://qgis.org/ubuntugis xenial InRelease
Obj:5 http://mirror.cedia.org.ec/trisquel.packages flidas-backports InRelease

Leyendo lista de paquetes... Hecho
Leyendo lista de paquetes... Hecho
Creando árbol de dependencias
Leyendo la información de estado... Hecho

No se pudieron instalar algunos paquetes. Esto puede significar que
usted pidió una situación imposible o, si está usando la distribución
inestable, que algunos paquetes necesarios aún no se han creado o se
han sacado de «Incoming».

La siguiente información puede ayudar a resolver la situación:
Los siguientes paquetes tienen dependencias incumplidas:
 python-qgis : Depende: libqgispython3.4.4 pero no va a instalarse
               Depende: libqgis-analysis3.4.4 pero no va a instalarse
               Depende: libqgis-core3.4.4 pero no va a instalarse
               Depende: libqgis-gui3.4.4 pero no va a instalarse
               Depende: libqgis-server3.4.4 pero no va a instalarse

 qgis : Depende: gdal-abi-2-2-2 pero no es instalable
        Depende: libgdal20 (>= 1.9.0) pero no es instalable
        Depende: libqgis-analysis3.4.4 pero no va a instalarse
        Depende: libqgis-app3.4.4 pero no va a instalarse
        Depende: libqgis-core3.4.4 pero no va a instalarse
        Depende: libqgis-gui3.4.4 pero no va a instalarse
        Depende: qgis-providers (= 1:3.4.4+24xenial-ubuntugis-1) pero no va a instalarse

 qgis-plugin-grass : Depende: qgis-provider-grass (= 1:3.4.4+24xenial-ubuntugis-1) pero no va a instalarse
                     Depende: libgdal20 (>= 1.8.0) pero no es instalable
                     Depende: libqgis-app3.4.4 pero no va a instalarse
                     Depende: libqgis-core3.4.4 pero no va a instalarse
                     Depende: libqgis-gui3.4.4 pero no va a instalarse
                     Depende: libqgisgrass7-3.4.4 pero no va a instalarse
                     Depende: grass760 pero no es instalable

E: No se pudieron corregir los problemas, usted ha retenido paquetes rotos.


texts: Question about quantum gis, and importing csv point data

texts: I'm trying to add some point data from a csv file (it's in decimal degree lat lon) and overlay  over road data that I already have added. The problem is that I can't get it to line up.  I'm sure I'm making some kind of mistake with the coordinate system, but I can't figure out what it is.  The data appears to the left and down miles away from my road data. When I import the same data in arcgis, I have to choose GCS_North_American_1983_HARN as my geographic coordinate system to make it line up, but that does not seem to be an option in quantum gis.  Any suggestions?


texts: How to create PRJ file for selfmade polygons in ArcMap?

texts: I have made a floodmap in ArcMap 10.4 and I would like to export it to ArcGIS Online. I tried, but it keeps saying ....shp can't be exported because it has no prj file. How do I create a prj file from normal shapefiles?


texts:
I am trying to calculate the population count (raster geotiff) lying under a polygon (vector layer). The appropriate QGIS library zonalstatisticsfb sum returns NULL. My QGIS python code is:
processing.run("qgis:zonalstatisticsfb", { 'COLUMN_PREFIX' : '_', 'INPUT' : 'C:/Users/Benedikt/test/isos.shp', 'INPUT_RASTER' : 'C:/Users/Benedikt/test/BY-AT.tif', 'OUTPUT' : 'C:/Users/Benedikt/test/out.shp', 'RASTER_BAND' : 1, 'STATISTICS' : [1] })

Input files can be obtained here in googledrive
The isos.shp attributes already contain the population count from openrouteservice isochrone population count (which I want to imitate). Why does QGIS calculate NULL? BTW: The QGIS Plugin creates the same result --> NULL.


texts:
I'm preparing a presentation of QGIS to release of the brazilian QGIS user guide. Thus, I would like insert some information about QGIS in global scale. Something like Tim's qgis-for-humanity presentation, but with updated information.
I'm very pleasure if you have something in this direction to send me.


texts: How to use the output of one run as input to the next in QGIS batch processing?

texts:
I'm trying to do a series of buffers using the processing batch interface, and I want to use the output of the first run as the input to the second (and so on). But I can only pick an existing file or open layer as the input, and since the last run hasn't finished yet it doesn't exist.

I'm pretty sure I can knock up a script to do this but it seems like the sort of thing you should be able to do using the batch processing interface so I'm wondering if I missing a trick?


texts: Problem loading specific WMS with leaflet

texts:
I try to add this WMS in a Leaflet map by using the L.tileLayer.wms.

https://glfc72.cfsnet.nfis.org/mapserver/cgi-bin/mapserver.cgi?map=/data/htdocs/glfc.cfsnet.nfis.org/mapserver/fire/public_html/ews_cmc_dev.map&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&CRS=EPSG:3857&BBOX=-20037507.842788246,-30240971.458386172,20037507.842788246,30240971.458386205&WIDTH=579&HEIGHT=874&LAYERS=fwi&STYLES=&EXCEPTIONS=XML&FORMAT=image/png&BGCOLOR=0xFEFFFF&TRANSPARENT=TRUE

And I get a 302 found error in Firebug. I am able to add the wms to ArcGIS Online without any problem but I would like to make it work in Leaflet because I do need a more personalized site. I am just starting with Leaflet so any kind of guideline about this would be very welcomed.


texts: Corrupted rasters after import into ArcGIS FGDB

texts:
I have been seeing a number of corrupted rasters recently after importing them into file geodatabases (FGDB).  The image below left shows a .tif prior to the import, while the image to the right shows the image after importing into the FGDB.  Could these errors be due to NoData handling?  What could be causing these errors?

Source TIFF info:

Columns and Rows: 5985, 7571
Bands: 1
Cellsize: 1,1 m
uncompressed Size: 21.61 MB
Source Type: continuous
Pixel Type: Unsigned integer
Pixel Depth: 4 Bit
Number of classes: 4
Storage Location: Google Drive


texts: spatialite : spatial index on geometries

texts:
I'm currently learning to use spatial access methods for query optimization. I'm going through the examples in spatialite cookbook and got stuck in here
http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/html/pp-adjacent.html

According to the tutorial, in order to find the pairs of populated places that their distance is < 1km we have to do this query

    SELECT pp1.name AS "PopulatedPlace #1",
      GeodesicLength(MakeLine(pp1.geometry, pp2.geometry)) AS "Distance (meters)",
      pp2.name AS "PopulatedPlace #2"
    FROM populated_places AS pp1,
    populated_places AS pp2
    WHERE GeodesicLength(
      MakeLine(pp1.geometry, pp2.geometry)) < 1000.0
      AND pp1.id <> pp2.id
    AND pp2.ROWID IN (
      SELECT pkid
      FROM idx_populated_places_geometry
      WHERE pkid MATCH RTreeDistWithin(
       ST_X(pp1.geometry),
       ST_Y(pp1.geometry), 0.02))
    ORDER BY 2;

Which throws an error because of no longer use of geocallback functions RTree***. I checked the update on this and this has been substituted by the use of virtual spatial index. We have to use a subquery of the form

    SELECT ROWID
    FROM SpatialIndex
    WHERE
       f_table_name = <table_name>
       AND search_frame = <some_geometry>

So i tested this

    SELECT pp1.name AS "PopulatedPlace #1",
      GeodesicLength(MakeLine(pp1.geometry, pp2.geometry)) AS "Distance (meters)",
      pp2.name AS "PopulatedPlace #2"
    FROM populated_places AS pp1,
    populated_places AS pp2
    WHERE GeodesicLength(
      MakeLine(pp1.geometry, pp2.geometry)) < 1000.0
      AND pp1.id <> pp2.id
    AND pp2.ROWID IN (
      SELECT ROWID
      FROM SpatialIndex
      WHERE
       f_table_name ='populated_places'
       AND search_frame = pp1.geometry)
    ORDER BY 2;

And the result set was awfully wrong

   PopulatedPlace #1    Distance (meters)   PopulatedPlace #2
              Ariano    0.000000             Ariano
   Campolongo Maggiore  0.000000             Campolongo Maggiore
   Campolongo Maggiore  0.000000             Campolongo Maggiore
   Campolongo Maggiore  0.000000             Campolongo Maggiore
   Campolongo Maggiore  0.000000             Campolongo Maggiore
   Campolongo Maggiore  0.000000             Campolongo Maggiore
   Campolongo Maggiore  0.000000             Campolongo Maggiore
              Ariano    0.000000             Ariano

Not giving me the pairs of places that their distance is <1km. Am I using a wrong spatial index subquery?  What i'm searching for is how to use properly the subquery to spatial index in order to obtain a number of rows of the table populated_places that satisfy a condition (pairs of places that are <1km away of each other in the above example)

NOTE: populated_places.geometry is of type POINT-2D and of srid=4236.

UPDATE: In case I do this :

    SELECT pp1.name AS "PopulatedPlace #1",
      GeodesicLength(MakeLine(pp1.geometry, pp2.geometry)) AS "Distance(meters)",
      pp2.name AS "PopulatedPlace #2"
    FROM populated_places AS pp1,
     populated_places AS pp2
    WHERE GeodesicLength(MakeLine(pp1.geometry, pp2.geometry)) < 10000.0
    AND pp1.id <> pp2.id
    AND pp2.id IN (234,235,236,237)

I get correct results but only for distances <10km. It seems as if there are not distances<1km which is false.


texts: ArcGIS cached map service sometimes pausing

texts:
I have published an ArcGIS service based on a MXD
I have generated the cached tiles on a shared NAS.
When I access the Map Service with the build in js admin viewer
(http://SERVER/arcgis/rest/service/....)

The map display quickly but, every now and then no tiles are coming for maybe 5 seconds, and then they all come quickly.

Why could that be?
Any tip / ideas welcome?
How to diagnose the problem?


texts: Trouble plotting Earth Engine images with Folium using Python

texts:
I am trying to plot an earth engine image using Python and the Folium package as described in the Python docs for Earth Engine. I tried running the example code given and I can't seem get the image to show up when saving as an HTML and opening in the browser. The getMapID() method always returns a token as empty string as in: mapID['token'] = ''. Does anyone know a way to get Earth Engine map tiles to show up on a Folium interactive map?

I also tried using the folium_gee functions by mccarthyryanc and ran into the same issues, the basemap shows up but not the tiles.

# Import libraries.
import ee
import folium

# Initialize the Earth Engine library.
ee.Initialize()

# Define a method for displaying Earth Engine image tiles to folium map.
def add_ee_layer(self, eeImageObject, visParams, name):
  mapID = ee.Image(eeImageObject).getMapId(visParams)
  folium.raster_layers.TileLayer(
    tiles = "https://earthengine.googleapis.com/map/"+mapID['mapid']+
      "/{z}/{x}/{y}?token="+mapID['token'],
    attr = "Map Data © Google Earth Engine",
    name = name,
    overlay = True,
    control = True
  ).add_to(self)

# Add EE drawing method to folium.
folium.Map.add_ee_layer = add_ee_layer

# Fetch an elevation model.
dem = ee.Image('USGS/SRTMGL1_003')

# Set visualization parameters.
visParams = {'min':0, 'max':3000, 'palette':['225ea8','41b6c4','a1dab4','ffffcc']}

# Create a folium map object.
myMap = folium.Map(location=[20, 0], zoom_start=3, height=500)

# Add the elevation model to the map object.
myMap.add_ee_layer(dem, visParams, 'DEM')

# Add a layer control panel to the map.
myMap.add_child(folium.LayerControl())

# Display the map.
myMap.save("myMap.html")


texts:
I have created a new field within one layer in Arcmap. and i have run some values and assigned each feature to a specific value inside the new field in the attribute table for the layer.

Now I am labeling for this layer and as i see, Arcgis have solved the overlapping issues for labeling the polylines features by its model.

I would like the Arcgis to solve the labeling overlapping issues for the polylines features by the weighted values, that i have created and added to the layer. Meaning that, i would like the Arcmap to give priority to each feature based on the value that it received from the field inside the attribute table.

I tried the Maplex, but i found out that you can choose the priority (High mod. low) for each layer and not for the features inside the layer.

Please can anyone give an idea how i can solve these problem. and if there is python or sql codes that can solve it, can you please give them to me.


texts: What professional mapping API would you use if your mobile app needed to seamlessly bridge both land and water with POIs?


texts: Including polygons as part of network for ArcGIS Network Analyst?

texts:
I am creating a pedestrian walking network of a downtown area.

I am digitizing line segments for paths and sidewalks, but I need to incorporate parking lots into my analysis.

So far, I have been creating multiple lines through parking lots to represent possible routes that a pedestrian might take through the lot.

This is tedious, and somewhat random.

I have my parking lots as a polygon feature class.

Is there any way for me to include these polygons in the network by snapping my walking network lines to the parking lot polygons?


texts: arcpy symbologyType is always OTHER, but I need UNIQUE_VALUES

texts:
I am trying to load a TIF file to the mxd using arcpy, then add a definition file to map colouring.

For one layer I have the problem, that the amount of values will always be different, sometimes it will be 0,1,2,3,4,5,6 sometimes it will be 0,7,8 and sometimes it will be 0,1,2,3,4,5,6,7,8,9,10,...,34

So i created in ArcMap an empty layer definition file: BEST_SERV_MO99.lyr
when added the TIF layer to the mxd, I am applying this. Now I would like to add all values except 0...
Also I would like to set the alpha to be 50%.

when I try the following code:

outRaster = os.path.join(path_now, tif_file)
addLayer = arcpy.mapping.Layer(outRaster)
self.log("addLayer: "+ addLayer.name)
arcpy.mapping.AddLayerToGroup(self.df, targetGroupLayer, addLayer, "BOTTOM")
layers = arcpy.mapping.ListLayers(self.mxd)
for layer in layers:
   if layer.name == tif_file:
       layer.name = tif2_file
       if os.path.isfile(symbologyLayer):
           arcpy.ApplySymbologyFromLayer_management(layer,symbologyLayer)
           self.log("layer: "+layer.name+" - symbologyType: "+layer.symbologyType)

I am always getting symbologyType: OTHER
The problem is, with OTHER, I can not add all values like this:

if layer.symbologyType == "UNIQUE_VALUES":
   layer.symbology.addAllValues()

also I have tried to do according this:

arcpy.mapping.UpdateLayer(df, updateLayer, sourceLayer, True)

but this did not have any effect, layer.symbologyType was still OTHER.

I am using ArcGIS for Desktop 10.3 with Python 2.7


texts: QGIS - Using Rule Based Labels With Related Tables

texts:
I have a two csv tables that have been related to eachother (RelatedResults):

Table 1 = Locations which is loaded as points
Which has the fields[ID] [Type] [X] [Y]

Example: Point_1, Frogs, LAT, LONG

Table 2 = Results which is loaded as just an attribute table
Which has the fields [ID] [Year] [ResultA] [ResultB]

Example:

Point_1, 2017, 2000, 4000

Point_1, 2018, 6000,150

I am trying to label the points based upon the year of sampling and every time I try to do it will only label using the result from 2017.

I have tried:

Filter

relation_aggregate( 'RelatedResults', 'count', CASE WHEN "Year" = '2018' THEN 'true' ELSE Null END )

Label Expression Generator

attribute(get_feature('Results', 'ID', ID), 'ResultA')

Not sure what I am doing incorrectly to get this to work as it will only ever show "2000" from the first year of results.

Any help would be amazing guys :). I have been stuck on this for a couple of days now. I need to produce concentration contours across a large areas for each year going back to 1985...


texts: RasterBrick from large rasterStack takes long to process in R

texts:
I have a list of 420 rasterLayer .tif files that I need to create a rasterBrick from.
I first create a stack:
library(raster)
rList <- list.files(pattern = ".tif")
rasterStack <- stack(rList)

Which is quick, but when I try make a brick from the stack, it takes really long.
rasterBrick <- brick(rasterStack)

I'm currently working on a cluster. Is this a job for parallel computing? I tried the following:
library(parallel)

numCores <- detectCores()
rList <- list.files(pattern = ".tif")
stacked <- stack(rList)
rasterList <- mclapply(1:length(stacked), mc.cores = numCores, function(x) brick(stacked[[x]]))

But I get the following error:
Error in seq.default(i, length(X), by = cores) :
  'by' argument is much too small
Calls: mclapply -> lapply -> FUN -> seq -> seq.default
Execution halted

I'm using all the cores available to me. Any suggestions please?


texts: Can I get the coordinates of country borders from Google Maps?

texts: I need to get the coordinates of a country's border (e.g., Iraq) from Google Maps. How can I do that?


texts: How to select multi-part polygons from a layer using PyQGIS

texts:
My script is -

# Find all multipart features in the active layer
l = iface.activeLayer()
iter = l.getFeatures()
geoms = []
for feature in iter:
    geom = feature.geometry()
    if geom.isMultipart():
        l.select(feature.id())
        geoms.append(geom)

print 'There are %i multipart features in this layer' % len(geoms)

I don't know the exact problem with this script.
Using this script in the Python console should select all multi part features in the active layer but in my case no multi section polygons fetched, while I can easily see multi section polygons exists in my polygon file.


texts:
I've written an open source Precision Agriculture application called AgOpenGPS that does everything from auto steer, to mapping, implement and section turn on/off, auto headland turns and more. It's used by farmers all over the world and everything is documented and free - a rarity in Ag. In ag we mainly use mapping in the area of a quarter section so lat/lon from the GPS is converted to UTM and the grid in meters is then used as basic cartesian coords for most everything.

One challenge I am having is if not farming right on the central meridian of course, the UTM grid declination at the zone borders causes a problem since driving straight north on the grid road should show 0, yet shows of course the error of a couple degrees. Heading can be determined by fixes of the grid, by VTG heading from GPS etc, however the two are only the same at the central meridian. If using a dual antenna for generating a stable heading, the auto steer will not work as the error is too great as the GPS heading vs the grid heading are a couple degrees off. Noticed the problem when AgOpen was being used alongside Trimble, Raven, Hemisphere etc. Years ago mfrs showed different headings, but since then they all aligned.

What I have done is calculated the Convergence angle using atan(sin(latitude)*tan(longitude- zone_north_meridian)). When opening or starting a new field, I take the integer of the northing and easting, store it, then any new fixes coming in are subtracted from those to end up with a manageable value to plot using openGL as a graphic display. (OpenGL does not like large numbers). So effectively this becomes my origin.

At this point i do a rotation matrix of every new incoming point based on the convergence angle and rotate that new fix and record its value as such. This does work, anywhere across the zone north is always north, going north on a grid road is north, the field edges line up with grid.

My question is,: Am I doing this correctly? To create field boundaries etc in Google Earth and import them, I do the same thing, rotate the imported points etc. Going back to google earth, I unwind them again. The major GPS manufacturers, systems costing $15 to 20k, do have it figured out, I am just wondering if I am going about it correctly.

ref: enter link description here


texts: How to plot Feature properties into a html div label using Leaflet and omnivore?

texts:
I have got some features in different GeoJSON files.

One of them includes Feature Collections (points and polyline) and the others just points.

I include those files into the Leaflet map using leaflet-omnivore plugin:

var track = omnivore.geojson("xxxx.geojson").addTo(map)

It is shown on the map as a polyline or a simple point, but I like to show a label including the "properties".

This is the .geojson file:

{
    "type": "Feature",
    "geometry": {
        "type": "Point",
        "coordinates": [-100.0947, 56.1581]
    },
    "properties": {
        "object": "aaaa",
        "psn": "lipsum",
        "lat": "lipsum",
        "long": "lipsum"
    }
}

I cannot find the properties into the variable track after the GeoJSON file is loaded.

I try to show a label with object,psn,lat, lon

All the examples I found show how to do that with a JSON object as a variable, but the goal is: how to identify the properties of the every feature once it is loaded through omnivore?


texts: Drawing with interior angle and distance

texts: In QGIS Azimuth and distance plugin draws with azimuth angle and distance. I want to draw with interior angle between consecutive lines. For the starting point and base station(reference point to plot 1st point) I have azimuth angles. Any plugin to do it?


texts: Syntax problem: Streaming GNSS Data from u-blox LEA-6H using str2str of RTKLIB not possible

texts:
Im trying to readout the data from my ublox LEA-6H recever over RTKLIB's  str2str command-line user interface (CUI). Im using rtkkib_2.4.2 on ubuntu 16.04.02 and the ublox receiver is connected over usb. im trying the following

sudo ./str2str -in serial://ttyUSB0:38400#ubx –out file://~/test.txt#ubx

i get the errormsg:

specify output stream(s)

how can i understand the documentation of rtklib. it says the following:


texts: How to specify TIME parameter in URL for WMS-T server in QGIS?

texts:
Does anyone know how to specify the time parameter for accessing data through a WMS-T server in QGIS?

I am trying to access fire data from FIRMS web server :(https://firms.modaps.eosdis.nasa.gov/web-services/#firms-wms-t)

The WMS-T GetCapabilities url is:
https://firms.modaps.eosdis.nasa.gov/wms-t/?SERVICE=WMS&VERSION=1.1.1&REQUEST=GETCAPABILITIES

How do I add the TIME parameter so that I can access the fire hotspots layers for specific dates?


texts: Selecting points using random field values in ArcGIS Desktop?

texts:
I have point layer, counties border area layer and I want to pick different number of points within different areas (counties) but from existing point layer - e.g. county A -> 10, county B -> 6, etc.

I'm wondering if code delivered by @Emil Brundage in Randomly subsetting X% of selected points using ArcPy? could be changed to use field values and area borders (like in Create Random Points tool) instead of one indicated count value to select random points from whole point layer. Because this code does only this.

Did I describe it clearly now?

def SelectRandomByCount (layer, count):
import random
layerCount = int (arcpy.GetCount_management (layer).getOutput (0))
if layerCount < count:
    print "input count is greater than layer count"
    return
oids = [oid for oid, in arcpy.da.SearchCursor (layer, "OID@")]
oidFldName = arcpy.Describe (layer).OIDFieldName
path = arcpy.Describe (layer).path
delimOidFld = arcpy.AddFieldDelimiters (path, oidFldName)
randOids = random.sample (oids, count)
oidsStr = ", ".join (map (str, randOids))
sql = "{0} IN ({1})".format (delimOidFld, oidsStr)
arcpy.SelectLayerByAttribute_management (layer, "", sql)


texts: How to export online content in QGIS print composer into static formats for offline use (eg: PNG, PDF)

texts:
I have a print composer layout which, amongst other things, has a HTML frame with the following code:

<class="hidden">URL</class>

<class="hidden">URL</class>

<class="hidden">URL</class>

This frame shows up fine on the screen after a short loading time (at least, most of the time, as from time to time it doesn't load and, hence, becomes empty), but my problem arises when I want to export all the atlas features using that frame into PNGs (or PDF, doesn't matter): everything is exported but that particular frame. I am assuming it may be due to a timeout or that it takes more time to load font awesome than the time QGIS takes in exporting, so my question is: Is there a way to prevent this from happening?


texts: How to split a single polyline shapefile into many shapefiles based on unique value using GeoPandas?

texts:
I have a single polyline shapefile with multiple features (row). I want to export each row to separate shapefile with using GeoPandas. But I got Geometry error in the code. This is my shapefile:

and this is the code:
Work_directory = r'P:/directory/Test/20200812'
folder_out = Work_directory + "/00_processing"
Shp_UniqueRun= Work_directory + "/SHP/UniqueRuns1.shp"
gdf = gpd.read_file(Shp_UniqueRun)

for index, row in gdf.iterrows():
    s = gdf.loc[index,'UniqueRun'] #use it for the name of output shapefiles
    shp_output = folder_out + "/" + str(s) + ".shp" #output shapefile
    NewSHP = gpd.GeoDataFrame(row)
    NewSHP.to_file(shp_output)

error:
AttributeError: No geometry data set yet (expected in column 'geometry'.)


texts: Get the count of no. of polygons after intersecting shapefiles

texts:
There are 2 shapefiles:
1. A polygon shapefile
2. Grid shapefile(which needed to be run in a loop)

Need to intersect the above 2 shapefiles and to get the count of no. of intersected polygons(which is grid wise).


texts: Loading Markers via XML and adding Infobox using Google Maps API V3?

texts:
I'm having a couple of issues that I can't seem to resolve so I'm going to throw this out there to the community.

I'm trying to create a Google Map that loads markers based off an existing XML file (I have that working). The issue I'm currently having is when the user clicks one of the Markers I would like the info window to show a "Category". Currently I can only get the Name, Lat, Lng to show up. When I try to show the category I get the info window but with no information.

The second issue I'm having is that it appears the Info Window is only working on one of my markers.

You can find an example of the map here:

http://www.officerage.com/maps/index.html

The example XML file that I'm using is located here:

http://www.officerage.com/maps/testxml.xml

Here is a copy of my code:

   <!DOCTYPE html>

    <meta charset="UTF-8">
    <title>A Basic Map</title>


texts: Plotting spatial data when two spatial objects have different CRS using R

texts:
I have a spatial polygon object and a spatial points object. The latter was created from xy latlong data vectors (called latitude and longitude, respectively) in a dataframe, while the former was simply read into R directly using rgdal. My code is as follows:

boros <- readOGR(dsn = ".", "nybb")
rats <- read.csv("nycrats_missing_latlong_removed_4.2.14.csv", header = TRUE)
coordinates(rats) <- ~longitude + latitude

At this point neither spatial object is projected. If I project these objects as follows:

proj4string(boros) <- CRS("+proj=lcc")
proj4string(rats) <- CRS("+proj=lcc")

Both objects are now projected, and both will successfully map with the plot() function as follows:

plot(boros)
plot(rats)

However when I try to plot them together:

plot(boros)
plot(rats, add = TRUE)

I get the first plot only, without the rats object superimposed over boros. However, and this is the big problem, I get NO error message, so I have been unable to determine where the disconnection is between these two spatial objects being able to speak to each other. Both commands run smoothly without error or warning, yet I am left with just the single plot. And when I check the projections of each object with proj4string() I get the same projection returned for each object.

I have now spent many, many hours over several days trying various ways of creating two spatial objects whose CRS and projections match such that they can be mapped together using plot(). Incidentally, one approach I took was to create a shapefile in ArcGIS for the rats object, which worked fine to create the file. But I am still left with the same inability of the two spatial objects to work together in R. I have tried many different projections for both objects, spTransform on both objects, nothing seems to work.

I have also included a dropbox link with the 2 data files I have described above: https://www.dropbox.com/sh/x0urdo6guprnw8y/tQdfzSZ384


texts: Using ArcPy Python Add-In Interval Caller to do auto-refresh?

texts:
I am using this code set called pythonaddin_interval_caller (written by Jason Scheirer) to do some map refresh every several seconds. The code was designed for arcgis python-addin. The benefit of this code is that it does the auto-check work without creating a separate thread.

    import tickextension.call_later
    import arcpy

    class Refresher(tickextension.call_later.CallQueue):
          def __init__(self):
              super(Refresher, self).__init__()
              self.call_later(self.refresh_view, 2)
          def refresh_view(Refresher):
              arcpy.RefreshActiveView()

          Refresher.call_later(Refresher.refresh_view, 2)

   refresh_ob = Refresher()
   refresh_ob.refresh_view()

This code is for python add-in. For some reason, I don't have permission to create add-in at this time. So I ran the script in arcmap python command window, I did two tests:

The original code only refreshed once, then, stopped. No error
I changed refresh_ob.refresh_view() to refresh_ob.refresh_view(refresh_ob), it throws the "takes exactly 1 argument (2 given) arcpy" error. I know its because I should not pass the object of self as the first argument. However, the program continued to run and did the auto-refresh every 2 seconds. The program works perfectly except this error message.

Are there any solutions to fix this?


texts: Concatenating names of the districts and the numbers of the sections in the polygonal layer of streets based on intersection

texts:
In QGIS 3.14 I want for a layer with street contours (polygons) to record the names of the areas through which the street passes and the numbers of cadastral plots through which it passes.
If a street intersects several objects, their names must be written in one field after a comma. How do I solve this problem?

if the area of the area falling into the boundaries of the street is less than the specified value, then its name is not taken into account


texts:
If you look at this site it has a main map window but it also has a small window map at the top right corner.

What is the technical term for the smaller window?
How is that window created?

I am currently using leaflet js and I wanted to insert that into my map here.


texts:
I am trying to get the x- and y-coordinate of a point layer in QGIS. So far, I have used the field calculator's function "$x" and "$y". I have read that easting coordinates should have 6 digits when using a UTM projection. However, when applying the functions mentioned, QGIS gives me a 7 digit easting coordinate. Does anyone know why this happens, and how I can obtain the 6 digit easting?
Projection: UTM 33 N
Example of coordinate: E: 1032624.32 N: 7862220.22


texts: Understanding the Heatmap tool values in QGIS

texts:
Context
I have been having a play around with the Heatmap (Kernel Density Estimation) tool in QGIS 3.4 and I realised that I do not really understand the output.
My understanding was that you choose the radius at which you want to analyse your point data with and the output will show you how many points can be found within a circle of chosen radius from any point in the output raster. I now know this to be wrong after looking at the results.
Example
Input
I have 3 points that I want to run the heatmap analysis (red circles) on and want to use a radius of 8 kilometers. I choose a pixel size of 100 (I understand the effect that varying pixels sizes can have on an output) and select the Quartic kernel shape (because it was default). Decay ratio is set to 0.

Output
The output can be seen in the below image. I have marked the raster values along a singular radius of each circle.

Question
What are these values telling me about that specific location in the heatmap and for a larger data set how could I make a heatmap where these values are meaningful to an audience who might not be as adept (I am clearly one of them!).


texts: Using ArcGIS API for Timeslider with Google Maps API?

texts:
Is it possible to use ArcGIS API for a Timeslider with Google maps API?

If this is not possible, then is there some other way I can use a timeslider outside of ArcGIS?


texts: How to assign correct flow direction to stream network polylines?

texts:
I would like to create a map showing the results of watershed analysis as well as the resulting stream network with flow direction.

I have a stream network which has been created from the flow accumulation raster. I also have a stream link raster which was created from the resulting stream network raster and the flow direction raster.

Once I convert the stream network or the stream links to vector polylines the resulting feature class does not retain the correct directions of the various segments.

What is the approach?

I have ArcGIS 10.2.1 Standard, Spatial Analyst, 3D Analyst

Red lines on the below example show the network I would like to show the map along with arrows indicating direction of flow:


texts: I found some maps from my state government that I wanted to I downloaded what I thought was a zipped up map perhaps as a bitmap. However, upon opening the zip file I found a large number of  .map .id .tab files. I take it that there is a program that can open all these files and create a map on my desktop. Can anybody tell me how to work with these files?


texts: Passing user-defined feature selection (by location) into geoprocessing widget in web app

texts:
I am creating a web based map with a custom script tool that is published as a geoprocessing service. I am using app builder to set up a geoprocessing widget, but am unsure of how to include a user selection.

I would like the user to be able to select features from a feature service layer in the web map (by drawing a polygon around the ones they want to select), and pass only the selected features into the geoprocessing widget that will run the script tool.

My initial thought is to have the user create a layer using a feature set and then write a select by location into the script to select features from another layer. However, I am unsure of what data type to use for the layer I want to make the selection on.

Is there an easier way to do this by using the select widget? Can you select features using the select widget and then pass only these selected features into the geoprocessing widget?


texts: Processing error when using PyQGIS in PyCharm

texts:
Running processing alghorithms works fine from the Python console in QGIS. I am using OSGEO4W (Python 3.7).

However when trying to run the same processing algorithm in PyCharm, I get the following error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:/OSGeo4W64/apps/qgis-ltr/python/plugins\processing\tools\general.py", line 105, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:/OSGeo4W64/apps/qgis-ltr/python/plugins\processing\core\Processing.py", line 183, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.

My preamble looks like the following:

import sys
import os

sys.path.append(r'C:/OSGeo4W64/apps/qgis-ltr/python/plugins')
sys.path.append(r'C:/OSGeo4W64/apps/qgis-ltr/python')
sys.path.append(r'C:/OSGeo4W64/apps/qgis-ltr')

os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = r'C:/OSGeo4W64/apps/Qt5/plugins'
os.environ['PATH'] += r';C:/OSGeo4W64/apps/qgis-ltr/bin;C:/OSGeo4W64/apps/Qt5/bin'
sys.path.extend([r'C:/OSGeo4W64/apps/qgis-ltr/python',r'C:/OSGeo4W64/apps/Python37/Lib/site-packages'])

from qgis.core import *
from qgis.analysis import *
from PyQt5.QtCore import *
from PyQt5.QtGui import *
from qgis.utils import *
from PyQt5.QtWidgets import *
import osgeo.gdal
import qgis.utils

QgsApplication.setPrefixPath(r'C:/OSGeo4W64/apps/qgis-ltr', True)
qgs = QgsApplication([], False)
qgs.initQgis()

from processing.core.Processing import Processing
Processing.initialize()
import processing
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
feedback = QgsProcessingFeedback()


texts: How parallelize the extract function for raster files in R?

texts:
I´m trying to parallelize a simple "extract" operation like this:

UseCores<-detectCores() -1
cl<- makeCluster(UseCores)
registerDoParallel(cl)
all.test.poly_period.test<-extract(stack.ts, test.poly)
stopCluster(cl)

stack.ts= stack of some hundred of images

test.poly= a shapefile with some polygons

But when I go to the task manager I can see that only one core is working.

I was wondering if anyone could give me a short advice about it?

The show() of the stack:

class       : RasterStack
dimensions  : 1596, 436, 695856, 186  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 284085, 297165, -2458245, -2410365  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=23 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : stack_MGM.1, stack_MGM.2, stack_MGM.3, stack_MGM.4, stack_MGM.5, stack_MGM.6, stack_MGM.7, stack_MGM.8, stack_MGM.9, stack_MGM.10, stack_MGM.11, stack_MGM.12, stack_MGM.13, stack_MGM.14, stack_MGM.15, ...
min values  : -0.26170975, -0.26808712, -0.27280286, -0.27959326, -0.28714523, -0.25714213, -0.22633411, -0.21442896, -0.20916982,  -0.20526922,  -0.17409700,  -0.16982903,  -0.17428745,  -0.18175986,  -0.18332276, ...
max values  :   0.9389138,   1.0646172,   0.8939760,   0.8914277,   0.8828803,   0.8622602,   0.8733163,   0.8952460,   0.9045494,    0.9184299,    1.0648371,    0.8974358,    0.8855312,    0.8626966,    0.8333160, ...


texts: Editing of vector layer programmatically with QGIS

texts:
I'm trying to edit vector layer using Python console and have couple of questions. Why this approach doesn't work and I have to use canvas.refresh? If I will use this code not in python console but in script file should I also use canvas.refresh?

# This approach doesn't work
layer = qgis.utils.iface.activeLayer()
layer.startEditing()
layer.moveVertex(7402946,6186167,0,0)

# This one works fine
canvas = qgis.utils.iface.mapCanvas()
layer = qgis.utils.iface.activeLayer()
layer.startEditing()
canvas.refresh()
layer.moveVertex(7402946,6186167,0,0)

What is the right way to perform such operation? I've found that I can make the same things using QgsVectorLayerEditUtils:

layer.startEditing()
eu = QgsVectorLayerEditUtils(layer)
canvas.refresh()
eu.moveVertex(7402946,6186167,0,0)


texts:
Default behavior of split in GIS software is mostly to update feature (row) for first part of split geometry, and insert copies for all other parts.
For many reasons this might not be desired behavior though, and instead delete original row, and insert all parts as new might be preferred.
In ArcGIS Pro - Catalogue - Feature Class Properties - Source tab there is option Split Model which can be switched from Update/Insert to Delete/Insert/Insert.
How to do this in QGIS?


texts: Invalid join selectivity error in postgis 2.1.3

texts:
I have 2 tables, one for points and one for polygons. I want to find points falling within theese polygons, so I quered:

select p.id, l.id
from points p inner join polys l on st_contains(l.geom, p.geom);

and i got invalid join selectivity error. I know that was issue in previous version of postgis but as far as I know it should be fixed in 2.1.3.
In my tables I have about 400.000 points and 1.000 polygons - maybe too many?


texts: Auto-recalculating M-values using ArcPy in ArcGIS Pro

texts:
I have M-enabled polylines, these lines are calibrated using different calibration points and some even start from zero. I want to drop existing M-values and then want to assign distance as M-value.
I can do this in ArcMap as you can see in the images below:
Dropping existing M-values.

Recalculating using Set as Distance.

I have many such polylines and I am trying to automate it using arcpy but could not find anything related to this.
Note: There is a tool  "Create Route", but it requires coordinate_priority which changes the direction of the line. I need to recalibrate line based on distance with direction preserved.


texts: Seeking GIS-related magazines and monthly papers

texts:
Can you recommend some magazines or other GIS-related papers (printed or digital) which are released on a monthly basis (or every 2-3 months);

I am looking to keep in touch with the GIS fast-growing world and I would prefer something other than books or websites.

Maybe some of you already have a subscription at some magazine.

I am quite interested in Open Source Software but not only.

This post would be the magazine equivalent to another one with useful information about books: What books, journals, and electronic resources are most valuable for expanding knowledge of GIS?


texts:
In QGIS 2.14 one can check if a rule gives results:

Is it possible to do this programmatically with python console? I would look for a way to iterate over all layers and there rules (nested rules!) and print those (layer_name:rule_name) which do not have results.


texts: Customizing legend with concatenated fields on QGIS

texts:
I prepared a map with a custom legend on QGIS 3.8.1, combining a classification of several concatenated fields. However, when creating the layout compouser and including the legend, the descriptive text of the legend did not keep the same line numbers.

What can be done to make legend able to keep the same number of lines defined in the layer?


texts: QGIS GPS position after entering attributes

texts:
I would like to be able to fill a form and, before, during or after, be able to cath a GPS point.
I already know the GPS extension but I have to catch the point and enter attributes after... Idea how I can do that?

Martin F


texts: Alternative cartographic software programs?

texts: I am looking for an alternative software product for cartography other than ArcGIS and Adobe Illustrator.  The graphic quality in ArcGIS outputs is substantially poor when compared to outputs in illustrator, yet illustrator is not a true cartographic software package.  Any open source?  Any extensions to ArcGIS or Illustrator that would enhance performance?  Read a couple referrals to MAPublisher in other forums, yet this package is too expensive ($2000+).


texts: Working a PyQGIS Script in cmd, outside of the QGIS python console

texts:
I have a perfectly fine working script that runs within the QGIS Python Console and I would like to run the script outside of the software, preferably in cmd as a .bat, so then I can use task scheduler to run the task every minute. I am using QGIS downloaded with OSGEO4W64. My .bat looks like this and runs, until it encounters a variable within QGIS e.g. QGSRasterLayer or QgsPrintLayout:
ECHO OFF
cd C:\Path\to\.py Script
C:\OSGeo4W64\apps\Python37\python.exe .\Full_Script.py
PAUSE

NameError: name 'QgsRasterLayer' is not defined

I know I have to do from qgis.core import QgsRasterLayer but every time I try to use this piece of code, it says that 'qgis' is not defined. I have looked into many forums on this site and have looked mainly into path and PYTHONPATH environment variables. My path looks like this:
C:\OSGeo4W\apps\Python37
C:\OSGeo4W\apps\Python37\lib
C:\OSGeo4W\apps\Python37\Scripts
C:\OSGeo4W\bin\qgis-bin.exe

And my PYTHONPATH looks like this:
c:\OSGeo4W64\bin\python

I have also tried using sys.path.append but I have had no luck so far. Can anyone help me out so I can sort out this problem and finally get my script working outside of QGIS, as I am struggling to find a way how to resolve this issue? Have I inputted the paths wrong or do I need to do something else, in order so that I can import qgis.core and be able to use all the elements (QgsRasterLayer etc) in my script, outside of QGIS.
EDIT: I currently have it working now but it comes up with the error: NoneType object has no attribute 'RasterLayer'


texts: Interpreting results of GWR - what do the values mean?

texts:
I used R to run a Geographically Weighted Regression (GWR) model.
The resulting model comes with a bunch of values that are given in a table. These are:

sum.w
X.Intercept, X.Intercept_se (standard error), X.Intercept_se_EDF
variable1, variable1_se, variable1_se_EDF
variable2, variable2_se, variable2_se_EDF,
variableX, ...
gwr.e
pred, pred_se, pred_se_1
localR2

Could you provide a short summary of what the calculated values mean?

Variable coefficient: I know that for each variable a coefficient is given, which shows the strength of the correlation. But are these values normalized and can be compared between variables? As an example: if variable1 has a coefficient of 5 in a certain area, does a coefficient for variable2 of 10 in the same area mean that its correlation is twice as high?

Standard error, EDF: What does the standard error tell? And what is the value EDF?

Sum, gwr.e, Intercept: I am also not sure what sum.w, gwr.e and the Intercept are for.

Prediction: The prediction values show what y-value is expected through the model.

Local R2: The local R2 shows the model performance at a local scale.


texts: ST_Intersect with ogr2ogr and Spatialite

texts:
I am trying to intersect Spatialite datasets with ogr2ogr, the OGR Virtual Format and SQL Syntax. I can't seem to locate whether the error is in my SQL query or my VRT file.

Data

rice_temperature.sqlite | Spatialite database containing 1 layer "rice_temp"

rice_elevation.sqlite | Spatialite database containing 1 layer "rice_elev"

rice_vectors.vrt

<OGRVRTDataSource>
<OGRVRTLayer name="temperature">
    <SrcDataSource>rice_temperature.sqlite</SrcDataSource>
</OGRVRTLayer>
<OGRVRTLayer name="elevation">
    <SrcDataSource>rice_elevation.sqlite</SrcDataSource>
</OGRVRTLayer>

All in the same directory I am also executing ogr2ogr from.

Query

ogr2ogr -f "ESRI Shapefile" intersection.shp rice_vectors.vrt -dialect sqlite -sql "SELECT ST_Intersection(rice_temperature.geometry,rice_elevation.geometry) AS geometry FROM rice_temp,rice_elev"

Errors

I tried a number of variations to the VRT file as well as to the SQL query itself. All errors are either no such table: rice temp or Cannot open datasource 'rice_elevation'

I have a feeling that I am close but can't wrap my head around what I am doing wrong - is it the SQL query or my use of OGRs Virtual Format?


texts: Python Error: “ValueError: need more than 1 value to unpack”

texts:
Can someone help me figure out what is not working in this piece of code:

import arcpy
import fileinput
import os
import string

from arcpy import env

env.workspace = "C:/EsriPress/Python/Data/Exercise08"
env.overwriteOutput = True
outpath = "C:/EsriPress/Python/Data/Exercise08"
newfc = "Results/newpolyline.shp"
arcpy.CreateFeatureclass_management(outpath, newfc, "Polyline")

infile = "C:/EsriPress/Python/Data/Exercise08/coordinates.txt"

cursor = arcpy.da.InsertCursor(newfc, ["SHAPE@"])
array = arcpy.Array()

for line in fileinput.input(infile):
    ID, X, Y = string.split(line, " ")
    array.add(arcpy.Point(X,Y))
cursor.insertRow([arcpy.Polyline(array)])
fileinput.close()
del cursor

I looked up some other questions/responses related to this error and they mention that perhaps the text file being read is not formatted correctly. Each line in my text file has an object ID, and an X & Y coordinates, with only one space in between each field (also, there is no extra line or space after the final object).

I'm wondering what else can be giving me this error.


texts: QGIS 3.2.1 - 3D rendering causes application to crash

texts: Whenever I try to render features on 3D map, QGIS crashes. Does anybody have any idea why? The 3D map works ok until I set 3D rendering for a layer.


texts: Converting mosaic dataset to TIN gives error?

texts:
I merged several DEMs together in on single DEM with mosaic dataset.
Now I wanted to convert this mosaic dataset in Vector with TIN. But when I use the tool (3D Analyst Tool) Raster to TIN comes an Error...

Have I to convert the mosaic Dataset first to another format? (How?)

Or do you have any other ideas to convert mosaic dataset to vector?


texts:
How to make a Heat Map of Distance between a set of parallel lines in QGIS?
This question is an extension of Question posted on the following link: Distance between multiple parallel lines?
I have a set of parallel lines (as in the image below), I am seeking to find the minimum distance of each line and create a heat map to visualize the distribution of distance between closest neighbors within extent of the layer in the direction of orientation of the lines.

Shapefile of the set of parallel lines is available from the following link: https://drive.google.com/drive/folders/1F8ZF6pZXB6RlqiASJ8m6J1f-LM9yEn4Z?usp=sharing


texts: ArcMap frequency table of field values that are above a certain value in another field

texts:
I am using ArcMap 10.5. I have a table of species list length from visits to sites over a period of years. I have a column for site, giving the grid reference of each site, and a column giving the number of species recorded on that visit. I want to create a frequency table that tells me how many visits where 2 or more species were recorded were made to each site.

Essentially I want to run the 'Frequency' tool on the 'site' field, but only calculate frequency for the rows where the 'species' field is >= 2.

How can I do this?


texts: Create Hexbin (honeycomb) grid using command line or Python

texts:
Qgis has this amazing Create Grid function. Is there an alternative for python / unix command line?

I've found Gdal_grid but it is not the same. The goal is to eventually create a hexbin heatmap.


texts: Possible ways to display maps offline, inside mobile hybrid app

texts:
I want to develop a mobile hybrid application to run on both Android & iPhone devices in which i want to display my own map data. I have the SHP files with the data i need. I have Geoserver and PostGIS for loading and serving the map data online. So far everything is fine.

My problem is that i still havent figure out the correct walkthrough for displaying the map data OFFLINE. Although i have read many answers here in this website, none of them seems to solve (as per my understanding my issue).

The goal is to display points of interrest on map. offline.
Some Points of interrest doesnt have address.Just Lon/Lat.
I want to cover these requirements:

How to offer the map data & the content offline to users device
How to offer to user routing options to navigate him to the destination from current location
How to add content data from a WordPress CMS, inside the app? (via JSON?GeoJSON?GeoRSS?)

Have been searching and looking many Tile services and many approaches to achieve the requirements, but still havent figure out the correct combination.

For the development of app I've been using PhoneGap XDK, Eclipse ADT. And for displaying the maps on mobile devices I'm using OpenLayers mobile.

Looking forward your suggestions or walkthroughs.

Thank you in advance,
with respect


texts: Spatially Joining Census Tracts to neighborhood shapefile?

texts:
I have a bunch of Census Tracts that I want to spatially join with larger neighborhood geographies so that the numeric fields of the Census Tracts are summed on the join.

A spatial join based on intersection won't work because some of the tracts intersect multiple neighborhoods. So I really want to join based on what neighborhood the centroid of the tract falls in.

However, when I run the spatial join tool with the center, it appears that center is only based on the target layer not the join layer. In older versions of ArcMap, I know I could run the Feature-to-Point tool to create a new feature of the tracts' centroids and then could run the spatial join that way. But this tool either doesn't exist anymore or isn't available with a basic license.

Does anyone have any workarounds? Or is there a better way to do this?


texts: Aggregate and split layer by location into multiple layer in ArcMap?

texts:
I'm pretty new to GIS, but hope I'm moving in the right direction.  I have a layer of the county map of LA, which is divided into communities and a layer of the building coverage of the county.  The building layer contains a large amount of data, so it takes a long time to render making visual analysis difficult.  It seems to me that the way to deal with this is to break the building layer into 353 much smaller layers based on community (there are 353 communities in the attribute table).  In that way I'm thinking that if a have a community_buildings layer for each community I can make them visible as needed without bogging down ArcMap.

It seems to me that what I would want to do is a spatial join on the communities layer to the buildings layer and then break that result into a 353 different layer based on communities, but I don't know how to accomplish the second part which is splitting one layer into multiple layers based on an attribute.

Below is the graphic of the data sets to get a feel for what I'm talking about.  Hopefully I'm moving in the right direction, but if I'm not can someone reorient me?


texts: Using requirejs with ArcGIS API for JavaScript v 4 API from CDN?

texts:
ArcGIS API for JavaScript Quick Start guides show "all-in-a-single-html" approach which is good for examples. However how about bigger apps?

I added comment on older issue (closed) at https://github.com/Esri/jsapi-resources/issues/32.

That works on 4.0, but I have requirejs import errors with the latest version today - 4.12:

Script error for "@dojo/framework/shim/Promise", needed by: esri/views/MapView
Script error for "@dojo/framework/shim/AbortController", needed by: esri/core/promiseUtils
Script error for "@dojo/framework/shim/array", needed by: esri/PopupTemplate
Script error for "@dojo/framework/shim/Set", needed by: esri/PopupTemplate
Script error for "@dojo/framework/shim/Map", needed by: esri/core/Evented
Script error for "@dojo/framework/shim/global", needed by: esri/request
Script error for "@dojo/framework/shim/WeakMap", needed by: esri/intl/date, esri/intl/number
Script error for "maquette", needed by: esri/views/overlay/ViewOverlay
Script error for "maquette-css-transitions", needed by: esri/widgets/support/widgetUtils
Script error for "maquette-jsx", needed by: esri/widgets/support/jsxFactory

I mean these require.config paths work on 4.0, but do not work on 4.12

paths: {
      'esri': 'https://js.arcgis.com/4.0/esri',
      'dojo': 'https://js.arcgis.com/4.0/dojo',
      'dojox': 'https://js.arcgis.com/4.0/dojox',
      'dijit': 'https://js.arcgis.com/4.0/dijit',
      'moment': 'https://cdnjs.cloudflare.com/ajax/libs/moment.js/2.24.0' // To resolve Script error for "moment/moment", needed by: esri/plugins/moment
    }

Repo with sample files for both version is at https://github.com/r-pankevicius/arcgis-js-api-with-requirejs


texts: In ArcObjects, how do I change the spatial reference of a Map (Data Frame) in an existing MXD?

texts:
I need to change the spatial reference of an existing MXD programmatically.  In ArcMap UI, I would take the following steps:  1) Open the data frame properties and go to the Coordinate System tab; 2) click on little globe button (Add Coordinate System) and choose New Projected Coordinate System; 3) choose Mercator; 4) enter custom name for the coordinate system; 5) change the Standard_Parallel_1 parameter to the appropriate value; and 6) Save/PressOK.  Ultimately, I need to automate this process over hundreds of MXDs.

I am really very confused about where even to start.  My research has turned up a lot of references to IProjectedCoordinateSystemEdit and related classes and interfaces and to reading/writing parameters to arrays with IParamter and related, but I am finding a lot of the information confusing.  So, which Interfaces/Classes should I start with, and what is the most concise way of accomplishing automation of the process I describe above?

My code currently opens the MXD to be processed in an IMapDocument and changes a suite of PageLayout and Map settings.  What is left is the changing of the spatial reference.  Thanks for any guidance you can offer.

UPDATE (@Hornbydd)  Here is the Python main body code.  (BTW, when accessing ArcObjects from Python, one thing that Python can't do is implement interfaces, so responding to events and extending ArcGIS isn't possible.  I have yet to find any com objects that it can't access, though.)

f = open(chartprops, 'r')

i = 0
for chrt in f:
    if i == 0:  #skip the header
        i = i + 1
        continue
    props = chrt.split(',')
    chartnm = props[0]
    scale = float(props[1])
    latorig = props[2]
    pixwidth = float(props[3])
    pixheight = float(props[4])
    if len(props) > 5:
        pprwidth = float(props[5])
        pprheight = float(props[6])
        pprorienttxt = props[7]
    else:
        pprwidth = pixwidth / 254.0
        pprheight = pixheight / 254.0
        if pprwidth > pprheight:
            pprorienttxt = 'LANDSCAPE'
        else:
            pprorienttxt = 'PORTRAIT'
    if pprorienttxt == 'PORTRAIT':
        pprorient = 1
    else:
        pprorient = 2

    #copy template chart MXD to chart output location
    #==== WARNING! ====
    # If the mxd already exists in the detination folder,
    # it will be overwritten
    #==================
    destmxd = outloc + os.path.sep + chartnm + '.mxd'
    shutil.copy(templt, destmxd)

    pMapDoc = NewObj(esriCarto.MapDocument, esriCarto.IMapDocument)
    pMapDoc.Open(destmxd)
    pPgLayout = pMapDoc.PageLayout
    pPgLayout3 = CType(pPgLayout, esriCarto.IPageLayout3)
    pPage = pPageLayout3.Page
    pPage.Units(esriInches)
    pPage.PutCustomSize(pprwidth, pprheight)
    pPage.Orientation(pprorient)
    pMap = pMapDoc.Map
    pMap.MapUnits(esriInches)
    pMap.SetPageSize(pprwidth, pprheight)
    pMap.ReferenceScale(1.0/scale)
    pMap.MapScale(1.0/scale)

    pMapDoc.Close()
    i = i + 1
f.close()


texts: ArcGIS-Calculate field : Give the first value to the others

texts:
I'm making a model builder to automatize some linear referencing projects but I have some problems with the field calculator tool.

My table contains 3 fields : ID | NOM | LIGNE and I'd like to complete the "LIGNE" field.

This row is numeric and it is composed by the 6 first number of the ID number 1.

I made a little picture to illustrate my problem :

(source: canardpc.com)
image

To resume, I'd like to complete my "LIGNE" field with the code "456789" which is the code line of the first ID.

Please... Could you help me ??


texts:
I have a problem with Orfeo Toolbox (OTB) function OrthoRectification.

I have a Pleiades product which is an extract (not a whole scene) of a small area in the western hemisphere south of the equator (so all negative coordinates). When I run the image through the OTB OrthoRectification I receive all kinds of problems with the output. When I order the app to produce output in cartographic coordinates it looks like the proper data, but is mirror flipped w.r.t. the east-west axis, so it's above (north of) the proper area AND has wrong resolution. See the image, which more or less depicts what happens:

When I attempt to use WGS84 (geographic coords) the geometry gets completely twisted and all I get is a weirdly stretched triangle.

To date I have tried all options of the Orthoreg function but no set of parameters can fix the problems. I have read the metadata that came with the image and the extent is OK there, so I assume it's some kind of mistake in OTB handling of coordinates, but can't be sure. Did you have a similar problem and if yes how did you overcome it?


texts:
I am conducting hydrological analysis for a catchment.
To make the DTM, I use mosaic to new raster tool.
I fill the DTM.
I make a flow direction layer.
I use the flow direction layer as input to the basin tool to show basins.

BUT... ArcMAP is merging two basins together when it should not.
See image.

I want to remove the basin to the west.

I tried reclassifying by increasing classes but the result is still the same with 32 (maximum) classes.

I tried this in MATLAB and it traced the catchment correctly. Maybe I am doing something wrong or it is simply an issue with ArcMAP not identifying basin divides properly? I have also tried three different DTM with different spatial extent and the result is the same as the attached image.

I am using ArcMAP version 10.4.


texts: Importing delimited text and applying style using PyQGIS?

texts:
I have managed to cobble together two Python snippets, one to import multiple delimited files, and another to apply a style. What I would like to do is combine these into one operation, and if possible assign a group for them (creating it ideally if it doesn't exist)....

The below is what I have so far:

import os.path, glob
from qgis.core import QgsProject, QgsVectorLayer
layers=[]
for file in glob.glob('C:/PathToFiles/*.QGIS'): # Change this base path and change backslashes to forward slashes
  uri = "file:///" + file + "?delimiter=%s&xField=%s&yField=%s&useHeader=no&crs=epsg:27039" % (" ","field_8","field_9") #(Delimeter, XField, YField) #1st column in file is column 1
  vlayer = QgsVectorLayer(uri, os.path.basename(file), "delimitedtext")
  vlayer.setFieldAlias(1,'X') #these order of columns go: 0, 1, 2 etc
  vlayer.setFieldAlias(2,'Y')
  layers.append(vlayer)

QgsProject.instance().addMapLayers(layers)

And:

from qgis.core import *
def applystyle_group(name):
    root = QgsProject.instance().layerTreeRoot()
    point = root.findGroup(name) #Find Group AP
    for child in point.children():
        if isinstance(child, QgsLayerTreeLayer):
            if child.layer().type()==0:
                child.layer().loadNamedStyle('C:/PathToStyle/POINT_CODE.qml')
applystyle_group("SPS Files")

As you can see, the group I am operating in is called 'SPS Files'.


texts: Maximum value of vector layer field in PyQGIS

texts: If I have a vector layer in QGIS, how do I use Python to find the maximum value which a given field has? I need to convert the values in one field into values between 0 and 1. I guess I therefore need to set it to val/maxVal. That's why I need to find the maximum value.


texts: Calculating MODIS 8-day composite GPP product as yearly

texts:
I want to use MODIS GPP 8-day composite product for an annual time series analysis.

How to calculate the annual GPP from the 8-day composite? Since the fractions represent the daily average rates of GPP for an 8-day period, do I multiply the sum of 8-day composite fractions of GPP by 8?


texts: What is Google source for United States topographic data?

texts:
I'm looking for as detailed as possible topographic data for a certain area and have found information displayed by Google Maps to be far more detailed than what is provided by 7.5 minute maps (from USGS).

Does anyone know the source of Google's topography data?

Specifically, I'm looking for data on the Los Angeles area.


texts: Can a QGIS GeoJSON vector layer be held in-memory?

texts:
When I add a vector layer to QGIS using data in GeoJSON format, data from the server are being displayed. The problem is that every change of the map extent (pan, zoomIn, ZomOut) is causing a request for the same data again.

Can the GeoJSON layer requested from the server be cached or held in-memory?


texts: QGIS Error : Could not load source layer for INPUT

texts:
When I try to clip a layer in QGIS 3.8.1 by using the toolbox clip function, it always returns the error Could not load source layer for INPUT XXX.tif layer not found.
When I manually drag and drop the layer to display it, it works.
I already changed the path to make sure it doesn't have something to do with the external HDD Im using or with letters (ä,ü,ö) in the path.
What else could be the reason?


texts: Trying to limit my shapemap to state boundary lines in R

texts:
I have a feeling I'm either trying to do this incorrectly, or I'm just not searching for the correct terms, but I cannot figure out how to do what I want to do.

I am trying to plot the ranges of fish species (shapefile #1) in Maryland/Delaware/Virginia superimposed over the state maps, then add major rivers (shapefile #2) in R.

I have no trouble getting everything plotted, but since both shapefiles are for the entire United States, they plot polygons/lines on my maps outside the MD/DE/VA region that I'm interested in.

Here's what I've done so far in R:

> library(maps)
> library(maptools)
> library(scales)
> mudsunfish=readShapePoly("~/Desktop/Maryland/Acantharchus_pomotis.shp")
> map('state', regions=c('maryland', 'virginia', 'delaware'), lwd=3)
> plot(mudsunfish, add=TRUE, col=alpha("skyblue", 0.3), border=FALSE)

This gives me the below image (not worried about the rivers in this instance):

Notice the blue shading in the upper right, where New Jersey and Pennsylvania would be. I don't want that. I can remove the shading in Photoshop, but with >120 species and rivers to add, I'd like to not have to do that.

How can I limit my shapefiles to the borders of the states? Should I be going about this with different packages? I'm sure someone else has run into this before, but I'm at a loss for what phrase(s) to Google to find my answer.


texts: I am the sole GIS manager for an organization of approximately 200 individuals that perform engineering assessments and collect other data to prepare for and respond to natural disasters. Most do not have any GIS experience. My supervisor wants to send them with an ArcMap template that they can fill in as needed. A topographic base map and style sheet are easy enough to send along, but I am wondering if it is possible to include a drop-down menu with preloaded options in the layout mode. The only way I can think of to do this is to use "Insert Object > Excel Document" and include this on the template.


texts:
I tried next command with geopackage, spatialite, geojson and others and they all process without any geometry. Shapefile just has DBF, geojson has , "geometry": null.

ogr2ogr -f "ESRI Shapefile" -skipfailures -progress -dialect sqlite -sql "SELECT NAME,TYPE_CODE AS TYPE,CLASS,UPPER_VAL || ' , '||UPPER_UOM || ' , '||UPPER_Code || ' , '||UPPER_Desc as UPPER,LOWER_VAL || ' , '||LOWER_UOM || ' , '||LOWER_Code || ' , '||LOWER_Desc as LOWER,LEVEL_CODE AS LEVEL,CITY,STATE,COUNTRY,CONT_AGENT as CON_AGENCY,COMM_NAME,SECTOR,ONSHORE,EXCLUSION,TIMESOFUSE|| ' , '||GMTOFFSET|| ' , '||DST_CODE as TIMES,REMARKS FROM Special_Use_Airspace_Pending" path\special_use_airspace.shp path\Special_Use_Airspace_Pending.shp

Source data is located at: http://ais-faa.opendata.arcgis.com/datasets/ac09bd0cf4164870bc7b68e60e8e09fb_0


texts:
I want to measure tree heights in the forest using drone photogrammetry.  So far it's easy to measure single trees height using fugroviewer in the pointcloud, but I would really like to measure whole groups of trees.

This picture below represents the tops. Is there a way to spot them out into the attribute table?

This is another image from before the treetops has been isolated.

I'm not a big GIS expert, so please be really specific in your explanation, if you have an idea.


texts: Filtering Image in Google Earth Engine based on metadata

texts:
I'm trying to create a raster layer of the USDA NASS Cropland dataset using only certain crop classes.  Is ee.Filter.inList only used for ImageCollections and not Images?  It seems more appropriate to import the layer as an Image, not an ImageCollection.

var cdl= ee.Image('USDA/NASS/CDL/2017').select('cropland')
var crops=ee.List([1,4,5,7,9,10, 31, 33, 36, 45, 63, 64, 86, 241])
var filtered=cdl.filter(ee.Filter.inList('cropland_class_names', crops))

This doesn't work because I get the error 'cdl.filter is not a function' on the last line. I need to export a new raster which only contains the crop types from a list.  I also tried filtering manually for each without success.

What is the proper way to filter my raster down to only these crop types?


texts: Extracting Lats and Lons from grib2 file

texts:
I am working on NCEP GFS data to plot certain parameters. I use python for plotting and wgrib2 for extracting data from grib2 files. There are no problem as such to extract parameter from grib2 data (and converting in to numpy arry). However, I need to get Lats and Lons also so as to plot using matplotlib library. The command line option to know grid (lats and lons) information does not seems to produce what is required i.e lat 90.0 to -90.0 with an interval of 0.5. I used this command to extract lats/lons data:-

wgrib2 -grid -d 1 ../GRIB_DATA/gfs_p50_grib_data/gfs.t00z.pgrb2.0p50.f000 -text -

and I received this data on terminal:-

-4.6
-4.6
-4.6
-4.6
-4.6
-4.5
-4.5
-4.5
-4.5
-4.5
-4.4
-4.4
-4.4
-4.4
-4.4
-4.3
-4.3
-4.3
-4.3
-4.2
-4.2
-4.2
-4.2
-4.1
-4.1
-4.1
-4.1
-4
-4
-4
-4
-3.9
-3.9
-3.9
-3.8
-3.8
-3.8
-3.7
-3.7
-3.7
-3.7
-3.6
-3.6
-3.6
-3.5
-3.5
-3.5
-3.4
-3.4

Can some one help me to guide in extracting lats and lons from grib2 data?


texts:
I created a heatmap containing population density data. Now I want to figure out the hotspots of this heat map. In other words: Where are the biggest population densities in my map?

My first thought:
Take a threshold and identify the biggest population densities. But somehow this method is a little too "random" for me, because its only me who decides where this threshold will be. I am searching for a method (e.g. automatic clustering) who will devide my density map into clusters and also gives an explanation why it took these clusters.

I am normally working with R, but I am also happy for solutions in another programm.

P.S.: For a better understanding I added the picture. White parts show high population density, black parts a low density. There are are some NULL Values.


texts: Ordnance Survey grid reference in QGIS print composer

texts: Is there a way of inserting Ordnance Survey grid reference in a print composer? I know that a grid can be created but I don't wish to insert a grid, I just to need to include the eastings and northings as a text string in the title block.


texts: Get EPSG code or CRS name from Proj4 WKT

texts:
I am trying to get the EPSG code from a proj4 WKT.

:“+proj=lcc +lat_1=46.183333333 +lat_2=47.483333333 +lat_0=45.666666667 +lon_0=-100.5 +x_0=600000 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs”

But I am unable to do so, I tried the same in QGIS but it also says user defined CRS. Although ArcMap is showig the CRS name correctly NAD_1983_StatePlane_Ohio_North_FIPS_3401_Feet. Is there any API or way to get EPSG code from a proj4 wkt.


texts: Getting historical administrative boundary GIS files for India?

texts:
I am currently working on a project for which I need to get GIS maps of India at a district-level in the 80's and 90's. India's district and state boundaries have changed considerably over the time period 1980-2011, and I need to map these changes over time.

Most of the versions that I found online are current (or based on the 2011 Census). Is it possible to find these older GIS files?


texts: How to get a Folium layer to display on top while not stopping lower layers highlight functionality

texts:
I'm trying to create a Folium map where the top layer is a boundary plot (fill_opacity: 0) of coarser structure to the Choropleth on the bottom level. I'd like the coarser level with no highlight functionally to show on top while the finer lower level to be shown on the bottom but it's highlight functionality to still permeate through and still work with the top layer not switched off with layer control. Is this posible?
My map:
m1 = folium.Map([54.5,-3],zoom_start=6.8,min_zoom=6,tiles='cartodbpositron',control_scale=True)

My top layer:
#makes boundaries plot
Boundaries = folium.GeoJson(
    cluster_stab_df,
    style_function = lambda x: {
        'color': 'black',
        'weight': 4,
        "opacity":1,
        'fillOpacity': 0,
    }).add_to(m1)

My bottom layer:
#sets style of Choropleth
style_function = lambda x: {"weight":0.4,
                            "color":'black',
                            "opacity":1,
                            "fillColor":colormap(x['properties']['gi']),
                            'fillOpacity':0.65}

#sets style of highlight
highlight_function = lambda x: {'fillColor': '#000000',
                                'color':'#000000',
                                'fillOpacity': 0.50}

#makes Choropleth map
Choropleth = folium.features.GeoJson(
    MSOA_stab_df,
    style_function=style_function,
    highlight_function=highlight_function,
    tooltip=folium.features.GeoJsonTooltip(
        fields=['msoa11nm','label','numneigh','gi'],
        aliases=['MSOA Name','Cluster Label','Number of Neighbors','Common Neighbors Proportion'],
        localize = True,
        style=("background-color: white; color: #333333; font-family: arial; font-size: 12px; padding: 10px;")
    )).add_to(m1)


texts: Extracting difference between two massive shapefiles using QGIS?

texts:
I have two massive shapefiles and I need to create a difference shapefile output between the two.  I have a layer where I have satellite coverage (orange) and then a layer where I have high resolution imagery (purple).  I'm looking to create a shapefile of only sections where I have satellite coverage but not high resolution imagery.  In my case these holes are where Landsat is used.

I'm not to used to doing anything on the command line but am happy to do so if it will work.

I believe the issue here is that the shapefiles are too big themselves (can contain >40k features).  None of the default difference tools work when using full shapefiles.  Will work for other queries though.

Here's a screenshot of what I'm talking about:

The final result that I need is the orange tiles with the areas covered by purple removed.


texts: How programatically activate data frame in ArcMap?

texts:
I want to:

create new data frame,
give it a name,
activate it and then fill with layers, etc. in ArcMap .NET custom command.

With this code I achieve 1. and 2., but don't know how to activate it:

// mxDocument is IMxDocument
IMap newMap = mxDocument.Maps.Create();
newMap.Name = "I can name it as I want: " + Guid.NewGuid();
mxDocument.Maps.Add(newMap);
//HOW TO Activate newMap???
mxDocument.UpdateContents();

New data frame is created, added to TOC but how to make it active?

mxDocument.ActiveView.FocusMap = newMap

gives System.NotImplementedException.

EDIT: IMaps.Create() +  IMaps.Add() is too primitive way to create a data frame, it won't work properly. See this link for a better description.


texts: Use ArcGIS Server REST geoprocessing service from the client without ArcGIS Web API's or ArcObjects

texts:
I'm checking to see if anyone knows how to execute and read the response from a custom Geoprocessor published on an ArcGIS Server as a REST web service without the use of ArcObjects or an ArcGIS web API on the client side.

I've built a custom geoprocessor service and successfully published it to ArcGIS server as a web service.  The service takes standard .Net objects as inputs and returns the same as outputs so there is no need for any ArcObjects on the client side.  I simply want to execute the REST web service by using standard .Net tools like System.Net.HttpWebRequest and System.Net.HttpWebResponse classes similar to what is shown here.

I've used the ArcGIS Server services with ArcGIS Web API's in the past, but doing it this way is new to me.


texts: Calculating relative position between geocoordinates and distances

texts:
I am working on a Unity application and am using the Haversine formula to calculate the distance between two points. I know the overall distance + the distance along the latitude and longitude separately given the formula. I know my current position and I know the heading (in degrees) using the compass. However, I am trying to conceptualize how to best get the relative positions of these data points.

Is the point in front of me or behind me as I rotate the camera and change the heading?

Is it to the left of me or right of me? (negative distance in world space)

I was originally planning to plot the points in World space based on their lat and lon. That still may be the best answer. However, the points are all pretty close together (~a few hundred meters) so the precision is pretty tight. I have been thinking of plotting the spaces along a 2-dimensional plane, offsetting them by their position relative to me and the heading location, hence my earlier question.


texts: Location of Workflow Manager Database References

texts:
Does anyone know where ArcCatalog references to Workflow Manager Databases are saved?

I am using Workflow Manager in a Citrix virtual desktop. Every time I restart the desktop after closing it, the reference is lost and needs to be re-created. I have not been able to locate any references in my AppData folders, where other ArcCatalog settings are stored.

I'm talking about the Workflow Manager Databases section of ArcCatalog's catalog tree not the Database Connections section.


texts: Common variable naming convention in ArcObjects

texts:
Not sure if this is really that important, but I'm wondering if there is a list of "common variable names" for different ArcMap objects...  for instance, sometimes you see:

IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument
IMxDocument mxDocument = ArcMap.Application.Document as IMxDocument
IMxDocument pMxDocument = ArcMap.Application.Document as IMxDocument

Then for IMap it could be:

IMap map = mxdoc.FocusMap
IMap pMap = mxdoc.FocusMap

Etc. etc. etc.    I'm just curious to know if there is a list of commonly used variable names (ie use mxDoc for IMxDocument,  use map for IMap,  etc.)


texts: Pyqgis QgsTemporalProperty visible features

texts:
I have a project with Temporal Controller, I have a QgsVectorLayer which 'Start and End Date/Time',I would like to know what features are visible  in every moment
I use Qgis 3.14, thanks


texts: Entering grouped data for single feature in ArcGIS Desktop attribute tables?

texts:
I have single objects, polygons of different zones. Each zone has been surveyed multiple times - ie. Zone A (one polygon) may have been surveyed in 2004, 2008, and 2012 and each year has several data fields (population, temperature, etc). In R I would structure my dataframe using "group_by" commands so that Zone A was the common thread between each survey line making for easy queries. How can I group these data in a similar way in ArcMap 10.3?

Ultimately, I want to be able to create a map symbolized by the number of surveys that occurred in each Zone, but retain all information from multiple surveys. I do not want a layer for each survey, just a single layer containing all Zones and survey data. I have dozens of zones and >150 surveys, so joins seem inefficient.

I've been reading to try to figure this out but think I lack the ArcGIS terminology to effectively search for what I'm looking for.


texts: Knowing acceleration, gyro and compass data how to calculate north direction

texts:
I have a small video camera with an attached sensor. I can get from the sensor his acceleration, gyro and compass data. Knowing this data I need to rotate the camera so that it will look strictly north.

I know how to rotate the camera, but how to calculate the right angle?

How can I do it? Or where can I read about it?
Or, if this is not a right place to ask such a question, what is the right place?

EDIT:

All coordinates are in an Earth-fixed inertial frame.The compass and accelerometer readings are unit vectors, and the gyroscope readings are in radians per second


texts: Calculating the fractal dimension of a line segment

texts:
I have a shp file with hundreds of features that represent clipped segments of coastline. Each segment is made up of many vertices (dozens, perhaps even hundreds). (FYI, the segments were created by using circles 1km in diameter to clip the coastline at random locations. The clipped segments of coastline are quite convoluted, I am looking to evaluate just how convoluted each polyline is.)

I need to calculate the fractal dimension of each individual segment (polyline). In order to be able to compare them.

I am working with ArcMap 10.1, I also just got QGIS but am not proficient with this program.

I know what the mathematical formula is to calculate fractal dimension for a single line segment, but I have far too many to calculate them individually, and I do not know how to write a script that could do it for me.

Using the box counting method:
Each line segment (or segment of coastline) is covered by a sequence of grids of descending sizes and then two values are recorded for each of the grids: the number of boxes that intersect with the line, N(L), and the side length of the box, L.

Note that N(L) ∝ L^(-D)

The fractal dimension is calculated as the slope of the linear regression best-ﬁt line of log–log data.

Therefore:
D = (logN)/[log(1/L)]


texts:
I'm fighting with pyqgis with a very simple issue.
I have a raster layer loaded in my TOC and I'd like to return with python its NoData value.

I looked in the API but (and sure it is my lack) I cannot find a way to do that.

So I get my raster as an object in the console:

rl = iface.activeLayer()

and I tried what I found in some websites (see the list at the bottom):

rl.getNoDataValue() #[1]
Traceback (most recent call last):
  File "<input>", line 1, in <module>
AttributeError: 'QgsRasterLayer' object has no attribute 'getNoDataValue'

rl.noDataValue() #[2] old API
Traceback (most recent call last):
  File "<input>", line 1, in <module>
AttributeError: 'QgsRasterLayer' object has no attribute 'noDataValue'

#[1]: http://qgis.sourceforge.net/qgis_api/html/classQgsRasterLayer.html#a10

#[2]: http://qgis.org/api/1.8/classQgsRasterLayer.html#a558e7f4bb4abe6b95c6882fbe4d2e19f

Look into the QgsRasterLayer class dir(rl) but I cannot find a way to return the NoData values.


texts: How do I clip Geojson FeatureCollection with another Geojson polygon

texts:
I have a GeoJSON file that contains one FeatureCollection with a bunch of features that are polygons. I also have the boundary of the continental USA in another GeoJSON file. I'd like to clip the first GeoJSON file, leaving only those features that are inside the bounds of the continental USA, preferably using ogr2ogr.

I've tried this:

ogr2ogr -f geojson -clipsrc us_boundary.geo.json out.json quad_5_usa.json

But I get a lot of these errors:

ERROR 1: TopologyException: found non-noded intersection between LINESTRING (-75.5474 35.7657, -75.4695 35.6083) and LINESTRING (-75.4806 35.561, -75.4917 35.666) at -75.490025631336778 35.649775563901855

and the output file is empty. Thoughts?


texts:
I have a surface model raster layer.  Is there any way to select only the light colored cells?  Or just the cells with a certain color characteristic such as color value?  This is a DSM of crop plots and soil.  I am using QGIS.

This is what I am trying to achieve by the use of zonal statistics, field calculator, and so on.


texts: Print Service legend display null for polygon features

texts: the print service of application works well, but  after converting the map area to pdf , the legend shows null filled values for polygon features(as in figure for Build up area,Forest dense area). I am using Map service to display those polygon features. but the polygon area that occupied for the category is not filled with relevant color value. I am trying to find a solution for this around two weeks still no luck. I am fed up of this. I use CMV esri javascript api framework to build my web application.


texts:
I'm new to QGIS and I'm trying to create a statewide map using specific regions, which break down by county. The steps I'm using are:

Vector Geometry: Fix geometries
Vector Geometry: Delete holes
Vector Geoprocessing tools: Dissolve by region

This works, but one county still has an outline on the map, even when the region map is the only layer (though if I check the attribute table, it is combined with its region).

In the attribute table, Shape_Leng, Shape_Le_1, and Shape_Area are all 0 for only this county, so I'm guessing this may be part of the issue. Is this likely the cause of the problem? And if so, how can I fix it? I got the shapefile from the SOS of the state and I don't have any other source.


texts: Estimating hydraulic head using ArcGIS Desktop?

texts:
I am trying to estimate hydraulic head using SRTM DEM. My objective is to calculate hydropower potential using the hydraulic head and inflow among others.

How do I get an estimate of head?

I am using ArcGIS Desktop 10.


texts: Geohash length to be reversible to LatLong using PostGIS

texts:
Decode/encode Geohashes are easy, and there are a "rule of the thumb" for estimat mininal Geohashe length for each LatLong precision. The question is other: precision for revesibility,

How many digits of Geohash for a symetric GeoURI with P decimal digits of precision? There are a (statistically) reliable formula for this kind of estimation?
There are other tool for PostGIS to calculate with better  precision the decode/encode of Gehashes?

Example: Geo URI    geo:-15.793889,-47.882778, with 6 decimal places, can be transformed into 6vjynkucf (9 base32 digits)  that seems good... But to back into the same LatLong we need (theorically) 12 digits.

SQL example with OpenStreetMap points

CREATE or replace VIEW vw_point_geohash_diffs AS
SELECT *, abs(orig_x-x) diff_x, abs(orig_y-y) diff_y
FROM (
SELECT orig_x, round(st_x(pt),6) x, orig_y, round(st_y(pt),6) y
    -- the 6 is "6 decimal digits at LatLong coordinates"
  FROM (
    SELECT osm_id, st_centroid(ST_GeomFromGeoHash( st_geohash(way,14)
                   -- here ...9, 10, 11, 12, 13, 14...
                 )) pt,
           round(st_x(way),6) orig_x, round(st_y(way),6) orig_y
    FROM planet_osm_point
 ) t1
) t2;
-- Precision checks:
SELECT n, sum_x, sum_y, round(avg_diffs/2.0,10) as avg_diffs, n_diffs_not0,
       round(100.0*n_diffs_not0/n,2)||'%' as perc,
       round(avg_diffs_not0/2.0,10) as avg_diffs_not0
FROM (
  SELECT count(*) n, sum(diff_x) sum_x, sum(diff_y)  sum_y,
       avg(diff_x+diff_y)  avg_diffs,
       count(*) FILTER (WHERE diff_x>0 OR diff_y>0) n_diffs_not0,
       avg(diff_x+diff_y) FILTER (WHERE diff_x>0 OR diff_y>0) avg_diffs_not0
  FROM vw_point_geohash_diffs
) t;

The theory is simple: merging all LatLong Digits in a integer number, 1579388947882778 and converting to base32, 1CSE8IJ6OOQ, results in a 11-digit number, that is the theoretic "mininal length". Converting each coordinate the length is 12:  base32(15793889)=F1VN1, 5 digits; and base32(47882778)=1DL8GQ, 6 digits; Imaging the worst case, base32(99999999)=2VBO7V, 6 digits. So 6+6=12.

Make sense, see data, a "S curve" for geohash_len X perc, with transition in len=12

geohash_len|  avg_diffs   | n_diffs_not0 |   perc  | avg_diffs_not0
-----------|--------------|--------------|---------|--------
10         | 0.0000020377 |      1786453 | 97.4%   |   0.0000020912
11         | 0.0000003716 |      1111317 | 60.6%   |   0.0000006130
12         | 0.0000000756 |       267843 | 14.6%   |   0.0000005173
13         | 0.0000000489 |       174916 | 9.5%    |   0.0000005127
14         | 0.0000000493 |       176324 | 9.6%    |   0.0000005128

Running all with n=1833326 points (OSM Brazil).

The main problem for "reversible precision"  seems the float arithmetic error at Geohash-14 or more. Same problem when we change round from  6 decimal places of LatLong to 5 places,

gh_len|  avg_diffs   | difs_not0 |  perc | avg_diffs_not0
------|--------------|-----------|-------|---------------
9     | 0.0000107338 | 1733044   | 94.53%| 0.0000113549
10    | 0.0000020370 | 679864    | 37.08%| 0.0000054930
11    | 0.0000003749 | 134963    | 7.36% | 0.0000050926
12    | 0.0000000836 | 30480     | 1.66% | 0.0000050300
13    | 0.0000000563 | 20502     | 1.12% | 0.0000050332
14    | 0.0000000570 | 20792     | 1.13% | 0.0000050293

The transition (S curve) is in len=10 or len=11, as expected.
 The problems arises after Geohash 13 digits, seems internal (floating point?).

PS: conservation of information  in translation or format convertion process have  its foundations in the notion of thermodynamically reversible process... The answer a/115501/7505  is not correct for this context of "reversibility".


texts: Why can't I use the output of Summary Statistics in my script?

texts:
I want to use summary statistics to dissolve duplicates in my table (no geometry) before performing a spatial join. The code works just fine and creates my output table in ArcMap but when I try to use it in my script I get an error that "GBAper1" does not exist. I have tried everything on my mind but can't find the problem.
What don't I see ?
  arcpy.Statistics_analysis(BldAreaPivot2, GBAper1, "afdak SUM;
  ingezonken SUM;loopbrug SUM;trap SUM;verdieping
  SUM;verheven SUM;zichtbare SUM", "CKEY")

  fields=arcpy.ListFields('GBAper1')
  for field in fields:
  print field.Object.name

When I write the same lines of code as specified in the tool help, I get a syntax error but the code above works as I get my outputs directly in Arcmap.
   GBAper1 = folderFoot + case + "_perGba1"
   arcpy.Statistics_analysis(BldAreaPivot2, GBAper1, [["afdak", "SUM"],
   ["ingezonken", "SUM"], ["loopbrug", "SUM"], ["trap", "SUM"],
   ["verdieping", "SUM"],["verheven", "SUM"], ["zichtbare", "SUM"]]", "CKEY")

   fields=arcpy.ListFields('GBAper1')
   for fieldObject in fields:
   print field.Object.name

Changing 'GBAper1' works for Listing the fields. I get the fields
(OBJECTID,CKEY,FREQUENCY, SUM_afdak,SUM_ingezonken,SUM_loopbrug, SUM_trap, SUM_verdieping, SUM_verheven, SUM_zichtbare). So GBAper1 exists and has a field "CKEY". But when I try to join it to a shapefile I get the message:

"ERROR 000728: Field CKEY does not exist within table Failed to
execute (AddJoin)"


texts:
Is there any way to create a styled map using a .json file from something like maputnik and export the map as a PDF or high resolution .png?

I currently have QGIS, maperitive, and tilemill, but I don't see any way to style these with the .json file.


texts: QGIS2web Ordnance Survey Map (raster layer) shifted against the QGIS project

texts:
I am wondering, why the QGIS2web can't accept the Ordnance Survey map location in the interactive version.
The Ordnance Survey map is set at EPSG 27700 British National Grid sor does a whole QGIS project and other layers too.
Other layers come the same as they're visible in QGIS, unlike the Ordnance Survey map canvas, which is the raster one. This one is shifted a bit and remains the loanword of EPSG 3857 Pseudo-Mercator coordinate system, want I can see when entering the code in layers.js.
       var lyr_Area6_1 = new ol.layer.Image({
                        opacity: 1,
                        title: "Area6",

                        source: new ol.source.ImageStatic({
                           url: "./layers/Area6_1.png",
                  attributions: ' ',
                            projection: 'EPSG:3857',
                            alwaysInRange: true,
                            imageExtent: [-23488.879733, 6814799.233467,
                -22713.936608, 6815622.095812]
                        })
                    });

      How can I fix it?

UPDATE:
In the QGIS2web plugin, there is an option: Match project CRS. When I tick it, then my map is much more shifted than previously:

I tried to use the coordinate capture and set these coordinates manually, but without result, unfortunately. Is there any way to fix it?
Again, both project as well as my OS raster map and the layers are in EPSG: 27700 British National Grid format.


texts: OpenGeo displays layer in OpenLayers, but not GeoExplorer?

texts:
My test server updated today I think it was a PHP update, and since then I can no longer view data in GeoExplorer, only OpenLayers. The GeoExplorer can not see my local server.

Any hints?

Update: Firebug says:

   TypeError: e is undefined
...&&!Ext.isEmpty(f.title)&&(e.set("title",f.title),a.setValue(e.get(d))):c.remove(...

 GeoExplorer.js (line 1535)

Another error from Geoexplorer:

Non-existing source 'local' referenced in layer config.

Further information from the GeoServer log:

org.geoserver.wms.capabilities.GetCapabilitiesTransformer$CapabilitiesTranslator.handleLayer(GetCapabilitiesTransformer.java:899)
    ... 87 more
2013-07-19 13:16:12,458 INFO [geoserver.wms] - Request: getServiceInfo

2013-07-19 13:16:12,459 ERROR [geoserver.ows] -
org.geoserver.platform.ServiceException: Could not determine geoserver request from http request org.geoserver.platform.AdvancedDispatchFilter$AdvancedDispatchHttpRequest@1965791


texts:
Trying to offset a line by 10 degrees, rotating around an endpoint.  I've read up on the older plugins and the new Rotate Feature tool, but am having trouble. I am using QGIS 2.8.9 Wien, reinstalled just a few weeks ago.

The Rotate Feature(s) tool allows me to rotate a line around a rotation point (which I know how to move), but I can't find any way to enter a specified rotation angle, or even to see what my angle would be.  I've read things that make me think it may have once been there, but can't find anything that says how to use it or if it can still work this way.
I installed CadTools and have tried using it, but keep getting errors:
An error has occurred while executing Python code:

Traceback (most recent call last):
File "C:/Users/bhughes/.qgis2/python/plugins\cadtools\tools\rotateobjecttool.py", line 75, in rotateObject geom = cadutils.rotate(self.feat.geometry(), self.pl,    angle • math.pi / 180)
File "C:/Users/bhughes/.qgisZ/python/plugins\cadtools\tools\cadutils.py", line 127, in rotate p1 = QgsPoint(i.x() - point.x(),  i.y() - point.y())
Attribute Error: 'NoneType' object has no attribute 'x'

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.8.9·Wien Wien, 26c3ece

It tells me that the CRS is undefined and that it's defaulting to WGS84, but ALL of my layers are defined as EPSG 102640, and so is my project CRS.  I've even got OTF turned off until the CadTools creates a CadLayer temporary layer in WGS84.
If I create a project in WGS84, the tool gives me the same error, but it will draw a line that I can save.  I got it to work once in 102640, but couldn't save the line; haven't been able to recreate it.
I don't want WGS84!

This is addressed on the hub under Feature Request #1787 and Feature Request #3472.  They both mention rotating by an angle; sounds like a window is supposed to just magically appear...

Am I the only one who can't make it work? Am I stuck in a WGS84 world?


texts: Changing default field Input Value Raster for Zonal Statistics as Table?

texts:
I have a raster in ArcGIS 10 with the following attribute table:

When I try to use the tool Zonal Statistics as Table over some polygons, it calculates the statistics of the fields 'Value' and 'Count'.
My goal is to keep the column 'Count', and give the values of the field 'Projection'to the field 'Value'.
Then it will be possible to use Zonal Statistics with the real values (922, 859, 929...) and their 'Count'.

If I use the tool 'Delete field' on 'Value', all the numbers move a column to the left.

If I use 'Reclassify', then I could manage to do it giving values one by one, but as there are thousands of values, and many more layers, this path is unfeasible for me...


texts: Automate “add junction flag-->trace downstream-->update field”

texts:
I have two point feature classes in a geometric network, say featureA and featureB. FeatureB is downstream of featureA. And there are about 10 featureB points downstream of 1 single featureA point. I need to update a field in featureB base on a field in featureA after trace downstream by selection.

I can do this manually by place “add junction flag” on individual point in featureA, and trace downstream by selection. Then update the field in selected featureB points(fieldB2) based on the value in featureA field(fieldA1).

The problem is:
I have 2000 points in featureA. I would have to repeat above process 2000 times.  Since it’s a repetitive process, is there any way I can automate this.

If FieldA1 is not Null: place a junction flag-->trace downstream by selection
Update selected fieldB1: FieldB1=FieldA2