texts: Dissolving/aggregating/combining points by proximity using ArcMap?

texts:
I'm teaching a GIS lab where we have groups of students go out and collect GPS coordinates of security phones on campus using the Survey 123 app.  Because of the way the it is set up, multiple groups are surveying the same phones and uploading them.  However, each phone's GPS gives a slightly different location.  So we end up with having multiple records of the same phone a few meters apart from each other.  If they were spatially coincident, its easy to use the dissolve function, but this doesn't work when they are not co-located.

Is there a relatively simple way we could have intro level GIS students go about "dissolving" these points?

For reference, we have about 25 phones they are collecting data for.  We could have them select the first record of each manually if necessary, but I'd rather try to show them a more justifiable approach.

If I were tackling this problem myself, I'd just use python to calculate the mean center of the multiple records for each phone.  But I don't really want to make them deal with trying to use a custom Python script or anything.


texts:
My employer provides me with a generous annual training budget and I'd like to use it on one of the following two Python courses.  As a GIScientist, I would like to attend the course which will give me the most robust skillset for both day-to-day tasks and larger customized interactive mapping applications.  In the future I would love to be able to develop custom web-map apps using all open-source (QGIS, PostGIS, etc), but am currently wedded to ESRI.  I do not have a programming background but have experience writing Python scripts mostly using Arcpy.

Learning Tree appears to be a more generic programming course with specific mentions of implementing classes and objects, accessing SQL databases, GUI development in Tkinter, and utilizing Django frameworks.

Enthought is geared towards scientists (my field) and engineers with specific mentions of Numpy, SciPy, matplotlib, and pandas for computation, building GUIs with Traits, and Chaco and Mayavi for visualizations.


texts: I'm trying to implement the Swedish language in Geonetwork 4.0.1-0. It went fine until it was time to make the schema-specific translations in the codelist.xml-file (\geonetwork\core-geonetwork\schemas\iso19139\src\main\plugin\iso19139\loc\swe). No matter how and what I try my translations won't take effect. Does anyone have any ideas for me?


texts:
I am using QGIS for mapping a series of buildings, I got the UTM coordinates from Google Maps and imported them as a delimited text layer.

However when I search for the UTM Zone that I need it doesn't appear I need 14Q and the only one I can get is 14N.


texts:
Using Q GIS 1.8.0 - I am trying to create several custom CRS definitions for the purpose of rectifying several split areas of a census shapefile.
In the "Custom Coordinate Reference System Definition" window I seem to be limited to creating no more than 6 custom CRS. If I advance the arrows to 6 of 6 and enter a new definition it overwrites USER 100006: CRS. This is repeated on any of the other CRS definitions (1 of 6, USER 100000, etc...).
Can I only create 6???


texts: Calculate sequential ID field using ArcGIS Field Calculator?

texts:
In this example, previously posted as Calculating sequential numbers into sorted table using ArcGIS Desktop? :

How do you construct a Python code that would give the Seq field sequencing such as:

1963.......1
1963.......1
1965.......2
1966.......3
1967.......4
1967.......4
1967.......4
1968.......5

Etc....   ???


texts: In ArcMap how do I print to an image (jpg or pdf) with a transparent background?

texts:
I am trying to print a contour shapefile to an image (jpg) with a transparent background instead of white.  What I want to do is add this into powerpoint as a layer over the top of a basemap so I can see both the basemap and the contours.

I realize that there are ways around this.  Example, I could just print two different maps; one with the contours and basemap and one with just the basemap.  I also know that I could also use the remove background in powerpoint, but we had issues with some of the maps not picking out the right areas.

However, my boss is adamant that there should be a way to do this exclusively in ArcMap.  Anybody know of a way.


texts: How can I use fortify() to create a filtered R data frame from a shapefile?

texts:
I am in the process of building a ggplot choropleth map of population in administrative areas in Wales. I have downloaded the Boundary-Line data from the Ordnance Survey and extracted what seems to be the right shapefile (community_ward_region.shp). Using R, I have got as far as reading in the shapefile.

require(maptools)
shape <- readShapePoly(wards)
str(shape)

Which gives me this promising output:

Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 1690 obs. of  4 variables:
  .. ..$ NAME      : Factor w/ 1507 levels "Abbey Cwmhir",..: 969 90 111 200 441 477 1455 249 255 305 ...
  .. ..$ DESCRIPTIO: Factor w/ 4 levels "COMMUNITY","COMMUNITY WARD",..: 2 2 1 1 2 2 2 2 1 1 ...
  .. ..$ COMMUNITY : Factor w/ 858 levels "Abbey Cwmhir",..: 67 67 81 128 152 152 152 152 157 190 ...
  .. ..$ FILE_NAME : Factor w/ 23 levels "ABERTAWE_-_SWANSEA",..: 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "data_types")= chr [1:4] "C" "C" "C" "C"
  ..@ polygons   :List of 1690
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 259009 188524
  .. .. .. .. .. .. ..@ area   : num 1923892
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:1629, 1:2] 259413 259420 259427 259427 259432 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 259009 188524
  .. .. .. ..@ ID       : chr "0"
  .. .. .. ..@ area     : num 1923892

Now if I do this:

bar <- fortify(shape, region = "NAME")

I get a nice data frame called bar that looks pretty much as I expected:

> str(bar)
'data.frame':   4744053 obs. of  7 variables:
 $ long : num  302962 302970 302974 303013 303015 ...
 $ lat  : num  280066 280076 280078 280097 280105 ...
 $ order: int  1 2 3 4 5 6 7 8 9 10 ...
 $ hole : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ piece: Factor w/ 29 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ group: Factor w/ 1762 levels "Abbey Cwmhir.1",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ id   : chr  "Abbey Cwmhir" "Abbey Cwmhir" "Abbey Cwmhir" "Abbey Cwmhir" ..

However, this is a large data frame and ggplot runs out of puff when trying to display it. In reality I only want to look at one area at a time. It looks as if the FILE_NAME factor in the shape object is what I want as it mostly corresponds to counties and the major conurbations.

> unique(shape@data$FILE_NAME)
 [1] ABERTAWE_-_SWANSEA
 [2] BLAENAU_GWENT_-_BLAENAU_GWENT
 [3] BRO_MORGANNWG_-_THE_VALE_OF_GLAMORGAN
 [4] CAERDYDD_-_CARDIFF
 [5] CAERFFILI_-_CAERPHILLY
 [6] CASNEWYDD_-_NEWPORT
 [7] CASTELL-NEDD_PORT_TALBOT_-_NEATH_PORT_TALBOT
 [8] CONWY_-_CONWY
 [9] GWYNEDD_-_GWYNEDD
[10] MERTHYR_TUDFUL_-_MERTHYR_TYDFIL
[11] PEN-Y-BONT_AR_OGWR_-_BRIDGEND
[12] POWYS_-_POWYS
[13] RHONDDA_CYNON_TAF_-_RHONDDA_CYNON_TAFF
[14] SIR BENFRO - PEMBROKESHIRE
[15] SIR_BENFRO_-_PEMBROKESHIRE
[16] SIR_CEREDIGION_-_CEREDIGION
[17] SIR_DDINBYCH_-_DENBIGHSHIRE
[18] SIR_FYNWY_-_MONMOUTHSHIRE
[19] SIR_GAERFYRDDIN_-_CARMARTHENSHIRE
[20] SIR_YNYS_MON_-_ISLE_OF_ANGLESEY
[21] SIR_Y_FFLINT_-_FLINTSHIRE
[22] TOR-FAEN_-_TORFAEN
[23] WRECSAM_-_WREXHAM
23 Levels: ABERTAWE_-_SWANSEA ... WRECSAM_-_WREXHAM

Q. How can I select only a subset of the data from the shape object I extract from the shapefile? For example, only the POWYS_-_POWYS parts? If I can somehow include the FILE_NAME in the data frame that is created with fortify then I could easily subset the bar data frame but I don't know how to do that. Or is there a way to use fortify to extract only parts of the object?


texts: Is there raster equivalent to vector Intersect?

texts:
The intersect tool that is used on vector data cannot be used on raster data.

What tool can I use for my two rasters that does the same as the Overlay > Intersect tool?


texts: Change detection using DEMs of difference and error quantification

texts:
Im currently working on a project in one of my seminars at university ("Advanced methods in natural hazard research").
We are modelling DEMs with Agisoft Photoscan using pictures from flume experiments (fluvial geomorphology) and are trying to detect changes between two different time steps.
Now i got stuck with the error correction of the DEM of difference (DoD):

"For the interpretation of the result, a very fundamental step is still missing: as our SfM approach has errors, we need to incorporate these into our change detection analysis. Luckily, PhotoScan has a function that creates a report on the model that we calculated, which also has information on the errors of our model (File  generate report). Save this report and extract the total error of your model. As both of the models that you calculated have errors, in a change detection analysis, we need to propagate these uncertainties. While there are many ways to achieve this, we will use the standard “Gaussian error propagation” for this. Propagating errors in a DoD analysis gives us the so called “level of detection” (LoD). Values in our DoD that are below the value of the LoD cannot be interpreted as significant changes, as we cannot differentiate them from our analytical error. The LoD using Gaussian error propagation is calculated using the formula (1). Within this analysis, we will also incorporate a degree of confidence, based on a student’s t-distribution:

= √(1)^2+(2)^2, (1); (square root after t!)
--> vertical errors are known but i dont understand what t is supposed to be.

Where t is the critical value of a given confidence interval in a student’s t-distribution, 1 is the vertical error of the model at time step 1, and 2 is the vertical error of the model at time step 2.
For the report (see below), create a figure showing the difference model you calculated, including the assessment of the LoD. Make sure that your colour scheme shows positive and negative changes in different colours and also the LoD is represented (by transparency or grey colours). Also calculate the magnitude of positive and negative changes and display them somehow, in a histogram, table, or some other form."

I understand that i have to delete every elevation that is <= LoD from the DEM of difference and create a new, corrected DEM of difference. Am I right?


texts: OpenGeo Suite Dashboard shows tomcat "localhost:8080" on browser instead of GeoServer?

texts:
This is the first time I'm working with GeoServer. For this reason I have installed openGeoSuite 4.7.1 and Java SE 7.

The issue is that when ever I try to open the application "dashboard" (path: C:/Program Files (x86)/Boundless/OpenGeo/Dashboard ) it shows me tomcat "localhost:8080" on browser instead of GeoServer.

OS is Windows 10.

Should I change the port number?


texts: Symbology effects on polygons in QGIS 3.4.10

texts:
I am wanting to represent polygons as dashed lines without filling the interior of the polygon. But there is the common edge between polygons, and so the visual effect of the lines vary. Depending on the zoom level (scale) there are places that appear as dashed lines, in other places that appear as continuous lines. How could I eliminate these visual effects in QGIS 3.4.10?

Figure 01

Figure 02

I tried the solutions for this question: How to remove line segments from a polygon for borders on a map in QGIS?

I tried the geometry generator option and also tried to create a virtual layer with the st_boundary option, but in both cases double lines are generated in the contacts of the polygons. I didn't find any command that generates an SQL without duplicate lines in the contacts.
I also tried using the "force right hand-rule orientation" option but I saw no result.


texts: qgis legend with calculated field symbol

texts:
I am making a map with a bunch of features plotted on it, each with a two-letter combination used as the symbol and a color. The letters are unique for each feature, drawn from a column containing the letter combination for each feature.
(That is, this is not a "label" per se, but a font symbol with a field being used as the input).

It works fine on the map itself, but when I want to output the legend, the symbols aren't being shown correctly, in stead qgis just shows a placeholder A on every single line in the legend.

Is there anything I can do to make this work?
I don't want to have to print the letters manually, because I have more than 2000 unique features (letter combination + color).


texts: Problem in adding raster layer from arcpy

texts:
I am using the following code to add the raster layer using arcpy in arcmap

def AddRasterLayer(path,layerName):
    layer=os.path.join(path,'new'+str(layerName)+'.lyr')
    mxd=arcpy.mapping.MapDocument("CURRENT")
    df=arcpy.mapping.ListDataFrames(mxd,'*')[0]
    arcpy.MakeRasterLayer_management(Raster(layerName),layer)
    newlayer=arcpy.mapping.Layer(layer)
    arcpy.mapping.AddLayer(df,newlayer,'TOP')
    new_name=arcpy.mapping.ListLayers(mxd,newlayer,df)
    new_name[0].name=str(layerName)
    del layer

After using this function raster layer is added to my current data frame. But when I am using any other tools it is showing as in the image

it is showing all the tree structure of the layer. Why is it so ? I need only the layers name to be displayed here.


texts: Why are results all black using gdalwarp via QGIS?

texts: I'm trying to convert a bmp image in BNG to WGS84 using gdalwarp (via QGIS). It seems to be running without errors but the output is all black. What am I doing wrong?


texts: How to speed up a query in a virtual layer of QGIS?

texts:
I'm trying to speed up a query in a virtual layer.
I noticed that the query runs at the same time for a partial extension and full extent of the layer shown on the map. So I thought if I filtered the query for the canvas extension I would gain an acceleration in performance.

I thought of something like this but don't know if is possible:

SELECT
      s.id       AS id,
      s.codigo   AS codigo,
      s.geometry AS geometry,
      SUM(d.e1)  AS sum1,
      SUM(d.e2)  AS sum2,
      SUM(d.e3)  AS sum3,
      SUM(d.e4)  AS sum4,
      SUM(d.e5)  AS sum5,
      DT_atu     AS DT_atu
 FROM "lyr1" AS s,
      "lyr2" AS d
WHERE d._search_frame_ = s.geometry
  AND ST_WITHIN(d.geometry, s.geometry)
  AND ST_WITHIN(s.geometry, iface.mapCanvas())
GROUP BY s.id;


texts: Counting intersections between two line layers using pyqgis

texts:
I want to count intersections between two line layers using pyqgis. First layer connects two points (Memorylayer_VL / Connections), the second one represents streets (lines_layer / Streets). I need to figure out how many times Connectionlines intersects with a Street. Note that the streets can have curves. So I also need to count, if a line crosses the same street more than once!
Thats where I am struggling to find a good logic how I could implement this. Anyone has a brilliant idea how this could be done? Thats what I have so far... It only counts if a street has been crossed, but not how many times.
    for tmp_line_feat in Memorylayer_VL.getFeatures():
        for streets in lines_layer.getFeatures():
            counter = 0
            if tmp_line_feat.geometry().intersects(streets.geometry()):
                counter = counter + 1
                attr = {5:counter}
                Memorylayer_PR.changeAttributeValues({ tmp_line_feat.id() : attr })

Example:

I need to know for each connection line, how many times a street has been crossed.


texts: Custom Geotool function that access a GeoServer's Datastore

texts:
I'm following this tutorial to implement a custom Geotools function:
http://docs.geotools.org/latest/userguide/tutorial/function.html

I want to implement a custom Geotools function that should access a GeoServer's DataStore. Could anyone please guide me to the code of connecting to a Datastore?

What I have tried is to access directly to the DB using the following code:

public static Connection createConnection() throws ClassNotFoundException, SQLException {
        Class.forName("org.postgresql.Driver");
        String url = "jdbc:postgresql://localhost:5432/spatialdbnew";
        Connection conn = DriverManager.getConnection(url, "postgres", "postgres"); //use postgres username to skip row-level security to enhance performance

        /*
         * Add the geometry types to the connection.
         */
        ((org.postgresql.PGConnection) conn).addDataType("geometry",
                (Class<? extends PGobject>) Class.forName("org.postgis.PGgeometry"));
        ((org.postgresql.PGConnection) conn).addDataType("box3d",
                (Class<? extends PGobject>) Class.forName("org.postgis.PGbox3d"));

        return conn;
    }

What I want is to access to the DB using the Geoserver's Datastore.


texts: How to place labels for varying projects?

texts:
I have a data set of features in a historic landscape. So if I make a map of it at a certain scale, I then need to move the feature labels around a bit—tweak them—to make it look nice.

But then if I do another  project to make another map of a smaller portion of the area—an enlargement at a larger scale—the labels again need adjusting—they need to be moved a bit to make it look good in that view.

However, if I then go back to my first project to revise that, all my label positions are messed up. Apparently the label location data goes with the original data. Seems like it would be better if label location adjustments like this were stored in a file with the project, so you could use the data for various projects at different scales.

Am I missing something? Is there a different way I should be doing this? The only resolution I can think of would be to make duplicate data sets for various projects using different scales.


texts:
I used worldclim 2.5 arc minutes data set for creation of ASCII file. According to worldclim 2.5 ac minutes represent spatial resolution of 4.5 km, therefore the cell size 4.5 x 4.5 = 20.25km2. but the properties of ASCII file showing cell size of 0.1 (x) and 0.1 (y), therefore the cell size is 0.01 square.

What is the correct one to calculate area of the desired pixel or cell in QGIS?


texts: How to draw a linear map programmatically?

texts:
I'm new to GIS programming and I'm using OpenLayers for my GIS project. I want to know how to draw a linear map programmatically between two locations from a real geographical map (may be OpenLayers). I think there should be a projected coordinate system for this. Is there any distance mapping algorithm to accomplish this task?

I do not mean adding a vector layer to an existing map. I want create a separate line map looking at a real geo map (Obviously a path connecting two locations).

Original map

Linear map


texts: Adding two rasters using map algebra in ArcPy?

texts:
I want to simply add two rasters within a python script using map algebra, which from what I understand is the equivalent to adding two rasters in the raster calculator.

If I do this using raster calculator and simply type "Raster1 + Raster2", it works perfectly. However, if I run this in a python script, I get error000732: Input Raster: Dataset Raster1 does not exist or is not supported.

Here is my code:

import arcpy
from arcpy import env
from arcpy.sa import *
raster1 = ("C:/raster1")
raster2 = ("C:/raster2")
raster3 = Raster("raster1") + Raster("raster2")

I set up the syntax based on the Esri help page for Map Algebra http://desktop.arcgis.com/en/arcmap/latest/extensions/spatial-analyst/map-algebra/an-overview-of-the-rules-for-map-algebra.htm

I'm not sure why my rasters are recognized within the raster calculator but not map algebra.


texts:
I have a set of rasters which I want to calculate mean of them,these rasters has some NoDate values,I have two approaches
first using this code in raster calculator I interpolate them

Con(IsNull("raster"), FocalStatistics("raster", NbrRectangle(5,5, "CELL"), "MEAN"), "raster")

Anyway this must work.But when I run it no changes happens.even by changing window size to 20 px nothing changes.

another approach is setting nodata values to 0 but It will affect on average results.
Another one I think I can do is ignore NoData values while avraging but for example when I have 4 rasters and one of them has NoData so in avrage it should devide into 3 not 4.

Anyway.Which one is better?any why I can not make first approach work?


texts: Is there a Geonetwork API to create a user

texts:
I would like to build a bit of custom code around geonetwork where I want to reuse existing database from gn. For my usecase it would be nice if users could register on a custom site and then use the credentials as their geonetwork login.

Now, instead of reinventing the wheel when writing all the neccessary tables like users, groups, email, etc. it would be helpful to have API access to those functionality.

The closest I found is this. So it's likely, something similar exists, but seems to be undocumented.

Does anyone know, how to access user creation functionality as is done when registering on the geonetwork-ui?


texts: ST_MakeSolid() creating an invalid solid from closed polyhedralsurfaceZ


texts: Adding Layers file to MXD using ArcPy- Code Error

texts:
In this code i try to add layer files- only if the first character in the name file, is a numbers to a MXD file, but an error happened:

import os,sys,fnmatch,shutil,arcpy

rootPath = r"G:\desktop\Project\lyr\taba\7650"
mxd = arcpy.mapping.MapDocument(r"G:\desktop\Project\Project.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
for root, dirs, filenames in os.walk(rootPath):
    for d in dirs:
        DIR = os.path.join(root,d)
        lyr_files = fnmatch.filter(os.listdir(DIR), '*.lyr')
        for i in lyr_files:
            if i[0].isdigit():
                addLayer = arcpy.mapping.Layer(i)
                arcpy.mapping.AddLayer(df, addLayer, "BOTTOM")
                mxd.save()

UPTADE:
I fix the code according to Hornbydd and still got the same error:

     for i in lyr_files:
            if i[0].isdigit():
                addLayer = arcpy.mapping.AddLayer(df, arcpy.mapping.Layer(os.path.join(root, i)))
                arcpy.mapping.AddLayer(df, addLayer, "BOTTOM")
                mxd.save()

this is the error:

Traceback (most recent call last):
  File "G:\desktop\abc.py", line 68, in <module>
addLayer = arcpy.mapping.Layer(i)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\arcobjects\mixins.py", line 389, in __init__
super(LayerMixin, self).__init__(lyrfile)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\arcobjects\_base.py", line 47, in __init__
for arg in args))
ValueError: Object: CreateObject Layer invalid data source
>>>


texts:
I have a shapefile of the worlds biomes

Terrestrial Ecoregional Boundaries Shapefile (TNC), 12/14/09, 50MB download from available from http://conserveonline.org/workspaces/ecoregional.shapefile

I would like to export it as a raster file; preferably a netCDF file with the MsTMIP format standard.

Ultimately, I could work with a netCDF that keeps metadata from the shapefile, or even a csv with no metadata. I can generate the .csv myself, but I expect that existing tools will do much of the work for me. Do they exist?


texts: Using OpenStreetMap Data in PostGIS database from ArcMap?

texts:
I've been trying for over a week now, and I'm not sure what I'm doing wrong and I haven't found the answer yet.

Target. I am trying to display and symbolize the OSM dataset over Jordan in ArcGIS for Desktop 10.3.

Software. Currently running PostgreSQL 9.3.5 and PostGIS 2.1 on a Windows 2008 R2 Server. ArcMap 10.3 is also installed on a Windows 2008 R2 Server. Along with ArcGIS for Server 10.3.

Steps after installation.

Use ArcGIS tool "Create Enterprise Geodatabase" and Name it jordan. (Works)
In PostgreSQL use CREATE EXTENSION postgis; (Works)
Use osm2pgsql using - osm2pgsql -s -c -H servername -d jordan -U sde -W -M -S .\default.style .\jordan-latest.osm.pbf to add the data to the PostGIS database (Works)
Connect to Database in ArcGIS (Works)
Register *_line database record with ArcGIS and Add OID (Works)
Add the feature in ArcMap (Works - Feature Geometry shows up)
Open Attribute table (Failed with error Could not load data from the data source. If you can correct the problem, press the refresh button to reload data. Possible problems can include bad network, invalid field, etc. Invalid SQL syntax [ERROR: syntax error at or near":" Line 1: select osm_id, access, addr:housename, addr:housenumber,... ::SQL state: 42602][jordan.sde.planet_osm_line])

I can not symbolize the data using the OpenStreetMap toolbox as well because of the same error.

I have tried changing the field names from ":" to "_" and the error goes away, but the next error states [ERROR: syntax error at or near","] and there is no comma in the field names so it must be inheriting it from somewhere else, and I have no clue where.

My ultimate goal of this project is to bring in the planet file from osm and symbolize it correctly.

What am I doing wrong, what other things can I try?


texts: Generating a map from a shapefile using ms4w mapserver

texts:
I am trying to generate a map from a shapefie. Through some initial searching, I wrote the code below but I am getting just a white page.

Below is my code:

 MAP
# To enter comments into a mapfile add a # symbol at the beginning of the line

NAME parcels1_map
STATUS ON
SIZE 800 800
EXTENT 489885 5450946 490904 5451637
UNITS METERS
SHAPEPATH "shapefiles"
IMAGECOLOR 255 255 255

WEB
TEMPLATE "template.html"
IMAGEPATH "C:/ms4w/Apache/htdocs/siya/output/"
IMAGEURL "siya/output/"
END

LAYER
NAME "parcel1_boundaries"
DATA states
STATUS DEFAULT
TYPE POLYGON
CLASS
NAME "parcels1_color"
COLOR 255 215 0
END
END
 # end of map file
END


texts: QgsApplication.processingRegistry() loads providers however is empty when called in within class function

texts:
I am in the process of building a standalone interface to QGIS to use the processing functions. In essence I am creating a mainWindow GUI with a Qgsmapcanvas. The GUI also includes some user inputs that feed into the algorithm parameters. So far I have had success in creating and running the GUI, map canvass and running the algorithms, however all separately. In the process of compiling the three I seem to be missing the processingRegistry files and as a result when I run an algorithm it can't find it in the registry (as it is empty).

The code snippet imports and loads everything, initializes the processing, runs the UI with the intention of pushing a 'run' button to execute the algorithm. My issue is it cannot find the algorithm. I have scattered 'printing the registry' throughout the code to see whats inside it at key points.

After the Processing.initialize() the full list of algorithms will be available but when it is queried within the class.execute_model() appears to be empty.
My initial thoughts and i think it revolves around my limited understanding of classes:

my execute_model function needs to be outside the class and call
the function from within the class?
I need to somehow pass the registry info into the class?
I've tried putting the Processing.Initialize() within the class
but that doesn't appear to work either.

My code:

import time
from qgis.core import *
from qgis.core import QgsApplication
from qgis.gui import *
from PyQt5.QtCore import *
from PyQt5.QtCore import pyqtSlot
from PyQt5.QtWidgets import *
from PyQt5.QtWidgets import QApplication, QDialog, QMainWindow
from PyQt5.uic import loadUi
from PyQt5.QtGui import *
from qgis.analysis import QgsNativeAlgorithms

# run QGIS Application
QgsApplication.setPrefixPath("C:/OSGeo4W64/apps/qgis", True)
app1 = QgsApplication([], True)
app1.initQgis()

# Run Process Application
import sys
sys.path.append("C:/OSGeo4W64/apps/qgis/python/plugins")  # Folder where Processing is located
import processing
from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
for alg in QgsApplication.processingRegistry().algorithms():
    print("{}:{} --> {}".format(alg.provider().name(), alg.name(), alg.displayName()))
from processing.tools import *

time.sleep(20)
# canvas_layers = []
# canvas = QgsMapCanvas()
# outputs = []

class IDC_Site_Selector(QMainWindow):

    def __init__(self):
        super(IDC_Site_Selector, self).__init__()
        loadUi("C:/IDCSS/IDCSSGUIv1.ui", self)
        self.setWindowTitle('IDC Site Selector v1')

        # check registry files
        for alg in QgsApplication.processingRegistry().algorithms():
            print("{}:{} --> {}".format(alg.provider().name(), alg.name(), alg.displayName()))
        #time.sleep(20)

        self.pushButtonRunModel.clicked.connect(self.execute_model)

    @pyqtSlot()
    def execute_model(self):
        results = {}

        # check registry files
        for alg in QgsApplication.processingRegistry().algorithms():
            print("{}:{} --> {}".format(alg.provider().name(), alg.name(), alg.displayName()))
        #time.sleep(20)

        # Load Master Files (to be replaced with a comboBox for LGA's in the future
        lots = ("C:/IDCSS/Model Data/Blacktown/Lot.shp")
        cadastre = QgsVectorLayer(lots, "vector", "ogr")

        # Add Lot Area Hectares
        alg_params = {
            'FIELD_LENGTH': 10,
            'FIELD_NAME': 'Lot_Area_Ha',
            'FIELD_PRECISION': 3,
            'FIELD_TYPE': 0,
            'FORMULA': ' $Area/10000',
            'INPUT': cadastre,
            'NEW_FIELD': True,
            'OUTPUT': "C:/IDCSS/Model Data/Blacktown/Lot_Area.shp"
        }
        feedback = QgsProcessingFeedback()
        results['LotArea'] = processing.run('qgis:fieldcalculator', alg_params, feedback=feedback)

app = QApplication(sys.argv)
widget = IDC_Site_Selector()
widget.show()
sys.exit(app.exec_())


texts:
I have a raw binary file containing color (RGB, 1byte per channel) with band interleaved by pixel (BIP) layout that I'd like to translate to a PNG. I could probably use some other non-GIS specific tool for this (recommendations welcome), but I'd prefer to use gdal_translate since it's already in my pipeline and has done a good job of translating other formats for me.

Is it possible to use raw files as input to gdal_translate?


texts: How to create projects in Tilemill Ubuntu Web service?

texts:
I have set up Tilemill on a Ubuntu server and successfully accessed it from browser using

TILEMILL_HOST:20009

I have already uploaded planet_osm data to postgresql/PostGIS DB in the same server machine. At first I have created a "Test" project where I have added a layer using the following information but the tiles were broken.

dbname=osm
host=<PG_DB_HOST>
port=<PG_DB_PORT>
user=<PG_USERNAME>
password=<PG_PASSWORD>

For the Subquery for testing pupose I used the following query

( SELECT way, name, amenity FROM planet_osm_point
WHERE amenity IS NOT NULL AND name IS NOT NULL
) AS data

Unfortunately there is no visualization in the web browser.

Basically I have created a custom design for an area in my local machine. I want to export MBTiles file to the server directly. I am using Tilemill for Windows locally where I have created that custom project. Now I want to transfer the MML project files and MSS Carto CSS files to the server and export MBTiles. The map changes will be frequent, this is why I need a solution where I can only upload the MSS and MML files to server and export MBTiles from Command line. The Local DB and Server DB Schema is similar. Now I have the following difficulties-

How to setup the project?
How to establish connection to PostGIS through Tilemill web service?
How to upload MML project file, MSS cart CSS files and other shapefiles and icons to server and which directory to upload?
How to export MBTiles? (Though I have read this Article. I am mentioning this question if there is any changes in exporting command while using Tilemill as web service)


texts: How can I get three different maps in one print composer?

texts: I'm using QGIS 2.6 and am struggling to get three different maps in one print composer. Better explained, one map with the county, the 2nd closer viewed with the town and the 3rd with the actual woodland. Its for a uni project however I can't seem to get beyond loading the 1st map. Iv saved it, I've clicked the lock layers but when it comes to reloading I just get a blank box reappear. Also, when I have the first map put on the composer, if I switch back to QGIS to load up the 2nd, print composer just closes. Any ideas?


texts: How to Select by Location and assign attributes using Python?

texts:
I have a polygon layer with a 'lot' field, and a points layer.

How can I use python to  select the points by location (to find what points are contained in what polygons) and assign the value in the 'lot' field from the polygon layer to the point that is contained in the relevant polygon feature?


texts: Is good database design less important for spatial databases?

texts:
I have a strong feeling that database design and normalisation often comes in second hand when dealing with spatial data.

With software costing a fortune and databases with over 100 fields tables I have to ask:

Is there good reasons for taking other considerations than normalisation when designing a spatial database?

I guess people will ask for examples, but that I can not give here, so my question is maybe more aimed for those who means that 100 fields is no problem and easier to maintain than a proper normalized design.

What are the arguments?


texts: Is it possible to have multiple subtype fields in a Feature class?

texts:
I am creating a feature class which will be used for data capture via a Geoform template in AGOL.  This has 17 fields.

I have set field 1 as a subtype with 3 options (3 types of conservation zones).
Field 2 contains the name of the conservation zones and I have set this up as 3 domains - one for each available option in field 1, each domain has around 30 different values.

I then also need field 5 to have its own set of subtypes (8 options) all of which will be linked to individual domains in field 6 (each domain holding around 30 different options)

Is this possible?

From the help pages on the ESRI website it appears that I am only able to assign one field to have subtypes, and any other "drop down menus" in the feature class  must be domains.  This means that everything is dictated by the subtype chosen in field 1.

The 2 fields I need as subtypes are completely independent of each other.

Does anybody know of a workaround for what I am trying to do?
Is there a better app than Geoforms in AGOL that will let me achieve what I am trying to do?

I am Using ArcMap 10.2.2 with an advanced licence and have an organisational AGOL account.


texts: Find every possible user to user travel time duration using pgr_dijkstra or pgr_kDijkstra

texts:
I have this below GetNearest function that will get every user's nearest node from ways (imported using osm2pgrouting).

CREATE or REPLACE FUNCTION GetNearest(NUMERIC, NUMERIC)
  RETURNS BIGINT AS
  $$
  SELECT id FROM ways_vertices_pgr AS vertices
  ORDER BY vertices.the_geom <-> ST_SetSrid(ST_MakePoint($1, $2), 4326) LIMIT 1;
  $$
LANGUAGE SQL;

Below is the query that runs and insert into another table the every possible user to user travel time duration.

INSERT INTO "UserTravelTime" ("FromLat", "FromLng", "ToLat", "ToLng", "Duration", "CreatedOn", "UpdatedOn")
SELECT C1."Lat" AS "FromLat", C1."Lng" AS "FromLng", C2."Lat" AS "ToLat", C2."Lng" "ToLng",
  (SELECT ROUND(SUM(cost)) AS Cost FROM pgr_dijkstra(
    'SELECT gid AS id,
         source,
         target,
         cost_s as cost
        FROM ways',
    (SELECT GetNearest(C1."Lng", C1."Lat")), (SELECT GetNearest(C2."Lng", C2."Lat")),
    directed := true)),
    NOW(),
    NOW()
  FROM
  "UserCoordinate" C1,
  "UserCoordinate" C2;

It works but it takes a lot of time to complete nearly 15 minutes, is there anything that I am doing wrong or could do better to make it faster?


texts: Gall-Peters or other cylindric equal-area projections in QGIS?

texts: Is it possible to use Gall-Peters or other cylindric equal-area projections like Hobo–Dyer in QGIS?


texts: ArcGis Engine app crashes when map contains class extension

texts:
I have an ArcGIS Engine 10.1 app. It has been running fine. It opens an .mxd on initialization. The .mxd contains 3 feature classes that have class extensions applied.

I just tried to make a change to the app for the first time in years and now I cannot get the app to load the .mxd.

The .mxd can be opened in ArcMap without a problem so I know the class extensions are installed fine on the machine. But when the ArcGIS Engine app loads it crashes without error or warning as soon as it loads the .mxd.

 axMapControl1.LoadMxFile("MyMXD")

If I remove the map layers from the .mxd that have the class extensions then it opens fine.

If I add the class extension library as a project reference and also have the class extension project open in visual studio I can step through the class extension code. What I see is the New sub being called and the module level variable being set, Once it gets to the end of the sub the ArcEngine app crashes, it doesn't call the New Sub on the next feature class in the map service.

Public Sub New()
    MyBase.New()
End Sub  'Crash after this line

I've just noticed an error being recorded in the output window aftwer the crash. When the app crashes after loading the first class extension the error

Unlicensed ARC session - terminating!

Is recorded in the Output window. The ArcEngine app successfully binds to a runtime and I successfully initialize a license with AoInitialize so not why this error is being thrown.


texts: Conversion problem of projection of MODIS LAI from sinusoidal to geographic

texts:
I am trying to change MODIS data's (LAI) sinusoidal projection to geographic. Actually I did it. But when I reproject it by ERDAS, the output LAI raster moved about 10 km to west side of actual shapefile and I tried it by ENVI, the output raster also moved about 20 km but to north side of shapefile. I am using the shapefile to check accuracy of the output raster. It was shape of lake.

Also I have tried it in R to convert projection using following script

library(sp)
library(rgdal)
library(raster)
Sys.setenv(TZ='GMT')

lai_real <-raster("E:\\LAI\\LAI_hdf\\MOD15A2_LAI_2014\\New Folder\\mod15a2.a2014169.h23_26.img")
projection(lai_real) <- CRS("+init=epsg:4326")
writeRaster(lai_real, filename="E:\\LAI\\LAI_hdf\\MOD15A2_LAI_2014\\New Folder\\mod15a2.a2014169.h23_26_geo.img", datatype='FLT8S', overwrite=TRUE)

But it doesn't work.

How can I convert the projection of the LAI from sinusoidal to geographic (lat/long)?

If my explanation is unclear, please visit the link that shared below. I hope then you understand what problem I have to convert sinusoidal to Geographic (WGS84)
https://www.dropbox.com/s/cucntnt8kh7m1dl/Hovsgol.PNG.


texts: NITF Metadata not displaying properly in OSGeo Shell. Works in FWTools

texts:
I'm using gdalinfo in the OSGeo Shell to display information on a directory of NITF's. For the majority of the NITF's, gdalinfo is saying the all the corner coordinates are 0, which is incorrect.

However, when running gdalinfo in FWTools Shell, all NITF's display coordinates properly.

This tells me it's not an issue with the metadata of the NITF's.

Has anyone noticed issues with support for NITF's in OSGeo Shell vs FWTools? It's seems like FWTools maybe reads header info differently then OSGeo Shell?


texts: Geometry errors when importing data from ESRI file geodatabase to PostGIS

texts:
When I run the following command to import by data from an ESRI fgdb to postgis...

ogr2ogr -f "PostgreSQL" -a_srs "EPSG:28355" -t_srs "EPSG:3857" PG:"dbname=tp_test user=postgres port=54321 password=xxxxxxx host=localhost" -overwrite P:\2012\183_TownPlanning_Symbology\Working\QueryLayer\TP_Query_Final_28355.gdb

...I get...!

  Warning 1: organizePolygons() received an unexpected geometry. Either
  a polygon with interior rings, or a polygon with less than 4 points,
  or a non-Polygon geometry.  Return arguments as a collection. Warning
  1: Geometry of polygon cannot be translated to Simple Geometry. All
  poly gons will be contained in a multipolygon.

I have checked geometry in ESRI and it reports no errors. I am running multipart to singlepart to see if this resolves the issues but this is crashing arcmap 10.1 as is dissolve. It ran a sort without issues and you can open the dataset.

I successfully imported the dataset earlier but then had to run the ESRI Dice tool to split the large polygons to no more than 15,000 vertices so that it processes faster.

I am new to postgis.

Is there a way to correct the issues in the new gdb or to split the polygons in postgis itself?


texts:
One of our ArcMap users is having a feature selection issue that I haven't run into before. I think the issues are all related, but when the user selects one feature, multiple features are being selected. I thought maybe it was the selection tolerance but that wasn't changed.

Then when looking into it a little deeper, I noticed when selecting in open space near the features, ArcMap is selecting the feature even though the selection tool is not near the actual features. As you can see in the image I am showing where I am clicking with the selection tool, and which features are being selected.

We are using ArcMap 10.3.1 and the features are stored in an ArcSDE SQL SQL Server 2008 Database .


texts: QGIS 2.18: Import Result Table of a query into QGIS as pointlayer

texts:
I have a query which outputs a table with an unique nodeID and x-,y-coordinates.

How can I import this query result table into qgis as an point-layer?

I have tried it with the DB-Manager (see below) and by WFS. But both hasn't worked out for me although I was able to specify a query.

The following shows you the SQL query I am running to get my result query table which I want to import in QGIS as point-layer:

SELECT "from" AS nodeID,
       st_x AS x,
       st_y AS y,
       priority AS type
FROM (SELECT "from",
             st_x(st_pointn(geom,1)),
             st_y(st_pointn(geom,1))
      FROM public.graph
      UNION
      SELECT "to",
             st_x(st_pointn(geom,st_npoints(geom))),
             st_y(st_pointn(geom,st_npoints(geom)))
      FROM public.graph
      ORDER BY "from") AS nodes
INNER JOIN junctions on nodes."from" = id


texts: Loading GeoJSON via AJAX after adding Layer to OpenLayers 3?

texts:
I'm trying to add a new GeoJSON layer to OpenLayers 3 that I have to fetch via AJAX. I want the user to click on a list of layers, have it add an empty vector layer, and then make an AJAX call and load the features into the layer. This is the function that's called when the user clicks on the layer:

function () {

    var styles = {
      'Point': [new ol.style.Style({
        image: new ol.style.Circle({
          radius: 10,
          fill: null,
          stroke: new ol.style.Stroke({
            color: 'magenta'
          }),
          cursor: "pointer"
        })
      })],
      'LineString': [new ol.style.Style({
        stroke: new ol.style.Stroke({
          color: 'green',
          width: 1
        })
      })],
      'MultiLineString': [new ol.style.Style({
        stroke: new ol.style.Stroke({
          color: 'green',
          width: 1
        })
      })],
      'MultiPoint': [new ol.style.Style({
        image: new ol.style.Circle({
          radius: 10,
          fill: null,
          stroke: new ol.style.Stroke({
            color: 'magenta'
          })
        })
      })],
      'MultiPolygon': [new ol.style.Style({
        stroke: new ol.style.Stroke({
          color: 'yellow',
          width: 1
        }),
        fill: new ol.style.Fill({
          color: 'rgba(255, 255, 0, 0.1)'
        })
      })],
      'Polygon': [new ol.style.Style({
        stroke: new ol.style.Stroke({
          color: 'blue',
          lineDash: [4],
          width: 3
        }),
        fill: new ol.style.Fill({
          color: 'rgba(0, 0, 255, 0.1)'
        })
      })],
      'GeometryCollection': [new ol.style.Style({
        stroke: new ol.style.Stroke({
          color: 'magenta',
          width: 2
        }),
        fill: new ol.style.Fill({
          color: 'magenta'
        }),
        image: new ol.style.Circle({
          radius: 10,
          fill: null,
          stroke: new ol.style.Stroke({
            color: 'magenta'
          })
        })
      })],
      'Circle': [new ol.style.Style({
        stroke: new ol.style.Stroke({
          color: 'red',
          width: 2
        }),
        fill: new ol.style.Fill({
          color: 'rgba(255,0,0,0.2)'
        })
      })]
    };

    var styleFunction = function(feature, resolution) {
      return styles[feature.getGeometry().getType()];
    };

    var layer = new ol.layer.Vector({
        source : new ol.source.GeoJSON({
            defaultProjection: "EPSG:4326",
            projection: "EPSG:3857"
        })
    });

    $.getJSON(global_config.services.getLayers + "?lname=fences", function (layerJSON) {

        if(layerJSON.features.length > 0)
        {
            var features = new ol.format.GeoJSON({
                dataProjection: "EPSG:4326",
                featureProjection: "EPSG:3857"
            }).readFeatures(layerJSON);

            layer.getSource().addFeatures(features);

            layer.setStyle(styleFunction);
        }
    });

    return layer;

}

The GeoJSON is in EPSG:4326 and the map is in EPSG:3857. The features are loading just fine into the layer, but they're not displaying and I can't figure out why. I'm thinking it has to do with the projection - though I have the proper projections set on all the features as well as the GeoJSON source layer. I've also tried having the default projection on the GeoJSON source layer set to EPSG:3857, because the features should be read in in that projection:

        var layer = new ol.layer.Vector({
            source : new ol.source.GeoJSON({
                defaultProjection: "EPSG:3857"
            })
        });

But still nothing is appearing.


texts:
I'm using a QGIS 2.8 in windows 10.
I have 2 scripts that are working separately but when I write both in just one script, it doesn't work. The 1st script is to load one raster to QGIs with a defined legend, the 2nd script is to compose a print layout and convert it to an image (png) in a output directory. My entire script is:

from PyQt4.QtGui import *
from PyQt4 import QtGui
from PyQt4.QtCore import *
from PyQt4.QtXml import *
from PyQt4.QtCore import QSize
from PyQt4 import QtCore
from qgis import core, gui

lyr = QgsRasterLayer("C:/PATRICIA/IG/TESTE/TEST/FWI_20150802.dat", "FWI_20150802")
s = QgsRasterShader()
c = QgsColorRampShader()
c.setColorRampType(QgsColorRampShader.INTERPOLATED)
i = []

i.append(QgsColorRampShader.ColorRampItem(8.5, QtGui.QColor('#2b9e27'), '8.5'))
i.append(QgsColorRampShader.ColorRampItem(17.2, QtGui.QColor('#fee717'), '17.2'))
i.append(QgsColorRampShader.ColorRampItem(24.6, QtGui.QColor('#ec7000'), '24.6'))
i.append(QgsColorRampShader.ColorRampItem(38.3, QtGui.QColor('#cc292c'), '38.3'))
i.append(QgsColorRampShader.ColorRampItem(100, QtGui.QColor('#75263b'), '100'))

c.setColorRampItemList(i)
s.setRasterShaderFunction(c)
ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1, s)
lyr.setRenderer(ps)
QgsMapLayerRegistry.instance().addMapLayer(lyr)

layers =iface.legendInterface().layers()
canvas=iface.mapCanvas()

for layer in layers:
    myFile = "C:/PATRICIA/IG/TESTE/Template.qpt"
    myTemplateFile = file(myFile, 'rt')
    myTemplateContent = myTemplateFile.read()
    myTemplateFile.close()
    myDocument = QDomDocument()
    myDocument.setContent(myTemplateContent, False)
    newcomp = iface.createNewComposer()
    newcomp.composition().loadFromTemplate(myDocument)
    newcomp.composition().refreshItems()

    c = iface.activeComposers()[0].composition()
    image = c.printPageAsRaster(0)
    image.save('C:/PATRICIA/IG/TESTE/OUTPUT/FWI_Arome.png','png')

My print layout composer appeared without layer! How can I solve this?

My current Script which print layout composer appeared without layer: `# coding: utf-8

#------------------------
from PyQt4.QtGui import *
from PyQt4 import QtGui
from PyQt4.QtCore import *
from PyQt4.QtXml import *
from PyQt4.QtCore import QSize
from PyQt4 import QtCore
from qgis import core, gui

lyr = QgsRasterLayer("C:/PATRICIA/IG/TESTE/TEST/FWI_20150802.dat", "FWI_20150802")
s = QgsRasterShader()
c = QgsColorRampShader()
c.setColorRampType(QgsColorRampShader.INTERPOLATED)
i = []

i.append(QgsColorRampShader.ColorRampItem(8.5, QtGui.QColor('#2b9e27'), '8.5'))
i.append(QgsColorRampShader.ColorRampItem(17.2, QtGui.QColor('#fee717'), '17.2'))
i.append(QgsColorRampShader.ColorRampItem(24.6, QtGui.QColor('#ec7000'), '24.6'))
i.append(QgsColorRampShader.ColorRampItem(38.3, QtGui.QColor('#cc292c'), '38.3'))
i.append(QgsColorRampShader.ColorRampItem(100, QtGui.QColor('#75263b'), '100'))

c.setColorRampItemList(i)
s.setRasterShaderFunction(c)
ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1, s)
lyr.setRenderer(ps)
QgsMapLayerRegistry.instance().addMapLayer(lyr)

layers =iface.legendInterface().layers()
canvas=iface.mapCanvas()

for layer in layers:
    myFile = "C:/PATRICIA/IG/TESTE/Template.qpt"
    myTemplateFile = file(myFile, 'rt')
    myTemplateContent = myTemplateFile.read()
    myTemplateFile.close()
    myDocument = QDomDocument()
    myDocument.setContent(myTemplateContent, False)
    newcomp = iface.createNewComposer()
    newcomp.composition().loadFromTemplate(myDocument)
    newcomp.composition().refreshItems()

    c = iface.activeComposers()[0].composition()
    image = c.printPageAsRaster(0)
    image.save('C:/PATRICIA/IG/TESTE/OUTPUT/FWI_Arome.png','png')`


texts: Displaying raster attributes on data driven pages

texts:
I have a feature class of 50 properties that are ranked. I have a photo for each property stored in a raster attribute column (in the feature class) called "Picture". I am trying to create a Dynamic text attribute for the picture column to use with my Data Driven Pages. I have struck out so far on getting the images to appear in the layout. Any ideas? Help would be greatly appreciated.

-Adam


texts: OpenStreetMap XAPI request for "Member Of"

texts:
I'm trying to find a way to request public transportation info from the OSM XAPI. Those nodes can be best described as:

highway = *

member of = ( route = bus )

Is that even possible to do over the XAPI or do i need to query it from the planet data?

this is what i tried but with no success

http://open.mapquestapi.com/xapi/api/0.6/relation[route=bus][bbox=100.63584242503583,13.659388026297398,100.65300856273045,13.66877069028344]


texts:
Not 100% sure if this is the correct place to ask my question so I apologise if it isn't.
Also, I have very limited experience with GIS systems, so please pardon my ignorance if some of my terminology is wrong, or if my question is actually very trivial.

I am writing an iOS application which will show a number of places on a map.
At my back end, I store the 'places' in the format:

Placename - string
Lat - decimal
Lon - decimal

I determine the current viewable area of the map (on the device) using the following Objective-C code:

MKMapRect mRect = self.mapView.visibleMapRect;
MKMapPoint neMapPoint = MKMapPointMake(MKMapRectGetMaxX(mRect), mRect.origin.y);
MKMapPoint swMapPoint = MKMapPointMake(mRect.origin.x, MKMapRectGetMaxY(mRect));
MKMapPoint nwMapPoint = MKMapPointMake(swMapPoint.x, neMapPoint.y);
MKMapPoint seMapPoint = MKMapPointMake(neMapPoint.x, swMapPoint.y);
CLLocationCoordinate2D neCoord = MKCoordinateForMapPoint(neMapPoint);
CLLocationCoordinate2D swCoord = MKCoordinateForMapPoint(swMapPoint);
CLLocationCoordinate2D nwCoord = MKCoordinateForMapPoint(nwMapPoint);
CLLocationCoordinate2D seCoord = MKCoordinateForMapPoint(seMapPoint);
NSLog(@"NW: %f, %f", nwCoord.latitude,nwCoord.longitude);
NSLog(@"NE: %f, %f", neCoord.latitude,neCoord.longitude);
NSLog(@"SE: %f, %f", seCoord.latitude,seCoord.longitude);
NSLog(@"SW: %f, %f", swCoord.latitude,swCoord.longitude);

This will give me an output of:

NW: 54.558793, -1.268541
NE: 54.558793, -1.254808
SW: 54.549659, -1.268541
SE: 54.549659, -1.254808

Its actually a box around my local bar :)

I want to pass the box coordinates to my backend system, which will return me information on any 'places' within the box.

What I require help with is the math to determine if a 'place' actually falls within the bounds of the box.

If anyone could give me an example of the math required, I would be very grateful.


texts: How to automatically run python code at the time when a mxd is saved?

texts:
I am trying to find a way to execute a few lines of (preferably python-)code whenever I save a specific mxd (my template.mxd)

I'm working with arcgis desktop 10.3.1 (advanced, all extentions)

What I am trying to achieve is the following:
I wrote a python script tool which enables the user to automatically export a lot of maps in pdf format. In order to do that, the user can choose certain basemaps from a value list. This list is generated in the validations of the script tool when the tool starts:

def updateParameters(self):
  mxd = arcpy.mapping.MapDocument(mxd_path)
  df = arcpy.mapping.ListDataFrames(mxd, "Detailframe",)[0]
  basemapList = arcpy.mapping.ListLayers(mxd,"basemap*",df)
  .
  .
  .

So now the user can choose each basemap which is integrated in the template mxd. But this slows up the tool startup (~20-30 seconds). Which bothers me.

Info: In the validations (def updateParameters) I made sure, that this list-generation runs only once instead of "whenever a parameter has been changed".

my idea is something like this:

The user alters the template-mxd by i.e. adding some new basemap layers and saves the mxd. At this point a second tool is triggered and lists all the layers with the suffix "basemap" and finally writes them to a simple textfile.
This text-file should be much quicker to access in the python script tool. So if the user opens the python script tool it reads out the previously generated textfile and gets all the latest layersnames from there. (the layernames is all I need...so a string or a list is fine)

I managed to do this in the label expressions in python code. But then the code of the second script runs every time when the user pans the map which massively slows up the template-mxd.

def FindLabel ( [text] ):

  mxd = arcpy.mapping.MapDocument(mxd_path)
  df = arcpy.mapping.ListDataFrames(mxd, "Detailframe",)[0]
  basemapList = arcpy.mapping.ListLayers(mxd,"basemap*",df)

  log = "c:/log/basemaptextfile.txt"
  txt = open(log, "w")
  for x in basemapList:
    txt.write('%s \n' % (x))
  txt.close()

  return

I read about arcgis addins which have similair functions which are called i.e. when a mxd is opened --> def openDocument(self):
How could I achieve that some code runs when the mxd is saved?

Of course this would ONCE slow down the template mxd...at the time, when the saving is done. But this would be much better than having to wait 30 seconds each time my tool starts.

Any ideas how to approach? I am not familiar with arcobjects, but I would give it a try if it isn't too complicated.

UPDATE:
the final construct shall have 4 components:

1) a template mxd which containes basemaps with the suffix "basemap". The user shall be able to alter this document by adding ne basemap layers or delete existing

2) python script tool which functions as a gui and reads out all the layers which have the suffix "basemap" and puts them into a value list, so the user can select which basemaps shall be exported to pdf. At this time this is achieved by listing the layers in the template mxd, but this takes too long. In the future this should be done by reading the layers from a txt.

3) a second code (maybe arcgis addin) which lists all the layers and writes them as a list or string into a simple txt

4) the python standalone .py which takes all these basemaps and creates pdf's from it(this works pretty well so far)

so what needs to be done is to find a solution for detecting if the user adds, removes or renames basemap layers in the template mxd and writing them to a txt. I thought this could be done by detecting if the template mxd was saved, because the user is supposed to save the template mxd before he runs the python script tool. But now I'm thinking this can maybe done with some functions from the arcgis addins. please feel free to suggest different approaches...


texts:
Is there a way to do a search for features that occupy more than one polygon in ArcMap?

I am working with a Parcel layer, and there are some cases where the parcels will take up multiple polygons (for example if they are bisected by a centerline).


texts: Newer version of QGIS doesn't display hillshade like in older version

texts:
I use QGIS 2.18 for looking at hillshades. I am trying to migrate to 3.18 but I am finding that the quality of the hillshades generated in the later version are not as detailed as in the older version. I can not find what I am missing. Could someone point me in the right direction to get this working right?
Side by side image comparison with 3.18 on left and 2.18 on right:


texts: Is there a way to automate field calculations for ArcGIS Server?

texts:
I have a field that I'd like to be concatenated with each new objectID when a feature is created in my feature service.

Workflow:

Feature added, already contains objectID.
Would like my Service Request field to be updated also on the DB side.
Service number = SRS + ObjectID = ex. SRS2012


texts: Earth Engine payload request exceeded on export

texts:
I am trying to export MOD13A2 data from Earth Engine. Export attempts have resulted in payload errors consistently, even after chipping away at the data volume. I tried removing clipping elements of the workflow based on this answer, but to no avail. It's weird because at ~24 million km^2, I don't understand why the 1km MOD product for only North America exceeds payload limits. I'm not performing complex image operations, and if it were an issue of raster size I'd expect Earth Engine to break the export into tiles. The only thing I can think of is that I flatten the imageCollection to a single image with (in this case) 20 bands - but I tried using ee.ImageCollection.toBands() and Tyler Erikson's approach here. These are both techniques I've used in other contexts with no problem.
The error,

Request payload size exceeds the limit:41943304 bytes

only occurs after the export has been initiated. That leads me to believe this isn't an image operation issue, since the first layer is drawn correctly in the plotting window. Surprisingly (to me) it seems like my Macbook goes into overdrive when I select "run" from the export dialogue box - not sure what my computer needs to be doing as this isn't a client-side operation as far as I know. But, the export dialogue freezes, the fan kicks on high and Google Chrome gets super slow; only after a few minutes does the error pop up and my computer returns to a normal operating situation.
Reproducible workflow below:
//// Data filters
// Create a region of interest
var ROI = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
  .filter(ee.Filter.eq('wld_rgn', "North America"));
Map.addLayer(ROI, {}, "ROI");
// Create a QA mask + clipping function
var masker = function(image){
  var mask = image.select('SummaryQA').eq(0);
  var maskedImage = image.updateMask(mask);
  return maskedImage;
};

//// Import data
// Import, filter by date, and apply QA mask
var MOD = ee.ImageCollection("MODIS/006/MOD13A2")
  .map(masker)
  .filterDate("2000-01-01", "2001-01-01")
  .select("NDVI");
// Look at the image collection
print("MOD",MOD);

// Plot the first image
Map.addLayer(MOD.first(), {min:0,max:10000}, "NDVI");

// Convert to single flattened image
var ndvi_Img = MOD.toBands();

// Export the data
// Export a cloud-optimized GeoTIFF.
// See https://developers.google.com/earth-engine/exporting
Export.image.toDrive({
  image: ndvi_Img,
  description: 'ndvi_Img',
  folder: "NoAm_MOD13A2",
  scale: 1000,
  region: ROI,
  fileFormat: 'GeoTIFF',
  maxPixels: 1.0E13,
  formatOptions: {
    cloudOptimized: true
  }
});


texts: Extract Census Blocks using Block Centroid that falls within Buffer in R

texts:
After creating buffer zones, I want to select a subset of Census Block IDs whose Centroids fall within the buffer zones. I have successfully built the data for all the Blocks and the designated Buffers, but I am confused by the literature (Geocomputation in R) on how to intersect the two datasets. I began with the following:

library(tigris)
library(sf)
library(tidyverse)
options(tigris_class = "sf")
options(tigris_use_cache = TRUE)
library(tidycensus)
census_api_key("123456")
library(haven)
library(foreign)
library(dplyr)
library(raster)

metro_blocks <- function(cbsa_code) {

  st <- states(cb = TRUE)
  cb <- core_based_statistical_areas(cb = TRUE)
  metro <- dplyr::filter(cb, grepl(cbsa_code, CBSAFP))

  stcodes <- st[metro,]$STATEFP

  tr <- rbind_tigris(
    map(stcodes, function(x) {
      blocks(x , year=2000)
    })
  )

  within <- st_within(tr, metro)

  within_lgl <- map_lgl(within, function(x) {
    if (length(x) == 1) {
      return(TRUE)
    } else {
      return(FALSE)
    }
  })

  output <- tr[within_lgl,]

  return(output)

}

# Pull dataset of blocks
blocks <- metro_blocks(47900)

blockst <- st_transform(blocks, 3488)

example <- data.frame(lat=numeric(), long=numeric())
example[1,] <- c(38.78309, -77.07458)
example[2,] <- c(39.44192, -77.40667)
example[3,] <- c(39.41428, -77.41073)

example_sf <- st_as_sf(example, coords=c("long", "lat"), crs=4326) %>% st_transform(3488)

# 1 mile = 1.6 km
# 5 mile = 8.0 km
example_5mi <- st_buffer(example_sf, dist=8000)

ggplot() + geom_sf(data=blocks[0]) + geom_sf(data=example_5mi[0], col="red", fill=NA)

Now I'd like to extract only the blocks (specifically column BLKIDFP00) that fall inside the Buffer Zones. Ultimately, I want a dataset that has every pair of BLKIDFP00, Buffer Points.


texts: 'Bring forward' Data Frame Extent Indicator in Arcmap

texts:
I want to 'order-> bring forward' the data frame extent indicator I have in an arcmap layout document. I want it to be on top of one graphic feature and under another.

This option is not available though, and other features in the document do not relate to it. How do I get it to appear 'on top'?!


texts: Where to find a version of libgeos-c1 compatible with libgeos-3.5.0?

texts:
This morning a number of Python scripts stopped functioning after a system upgrade. libgdal is no longer available, resulting in the following exception:

ImportError: libgdal.so.1: cannot open shared object file: No such file or directory

Trying to re-install the library by hand I am getting unmet dependencies regarding the libgeos-c1 package:

$ sudo apt install libgdal1h
Reading package lists... Done
Building dependency tree
Reading state information... Done
Some packages could not be installed. This may mean that you have
requested an impossible situation or if you are using the unstable
distribution that some required packages have not yet been created
or been moved out of Incoming.
The following information may help to resolve the situation:

The following packages have unmet dependencies.
 libgdal1h : Depends: libgeos-c1 (>= 3.4.2) but it is not going to be installed
E: Unable to correct problems, you have held broken packages.

In fact libgeos-c1 can be installed, but it would force the downgrade of libgeos, which would remove pretty much all other GIS software from the system:

$ sudo apt install libgeos-c1
[sudo] password for desouslu:
Reading package lists... Done
Building dependency tree
Reading state information... Done
The following extra packages will be installed:
  libgeos-3.4.2
The following packages will be REMOVED
  gdal-bin grass-core libgdal20 libgeos-c1v5 liblas-c3 liblas3 liblwgeom-2.2-5
  liblwgeom-dev libopenscenegraph100v5 libosgearth3 libosgearthannotation3
  libosgearthfeatures3 libosgearthqt3 libosgearthsymbology3 libosgearthutil3
  libqgis-analysis2.14.3 libqgis-app2.14.3 libqgis-core2.14.3
  libqgis-customwidgets libqgis-gui2.14.3 libqgis-networkanalysis2.14.3
  libqgis-server2.14.3 libqgisgrass7-2.14.3 libqgispython2.14.3 libsfcgal1
  libspatialite5 openscenegraph-plugin-osgearth postgis
  postgresql-9.3-postgis-2.2 python-gdal python-pyspatialite python-qgis
  python-qgis-common python-shapely python3-gdal qgis qgis-plugin-globe
  qgis-plugin-grass qgis-provider-grass qgis-providers
The following NEW packages will be installed
  libgeos-3.4.2 libgeos-c1
0 to upgrade, 2 to newly install, 40 to remove and 0 not to upgrade.
Need to get 0 B/454 kB of archives.
After this operation, 218 MB disk space will be freed.
Do you want to continue? [Y/n]

I have libgeos-3.5.0 installed from the UbuntuGIS PPA:

$ apt-cache policy libgeos-3.5.0
libgeos-3.5.0:
  Installed: 3.5.0-1~trusty1
  Candidate: 3.5.0-1~trusty1
  Version table:
 *** 3.5.0-1~trusty1 0
        500 http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu/ trusty/main amd64 Packages
        100 /var/lib/dpkg/status

If I understand correctly, I need to find a version of the libgeos-c1 package that is compatible with libgeos-3.5.0. Would anyone know where I could get it?


texts: Change FeatureCollection into Multipolygon

texts:
I have code that works with features that are multipolygons.
The problem is that I want to work with shapefile that I upload, but it's always FeatureCollection and my code fail for this.I have see nthere is a code to make geometry into featurecollection but I want to other direction.

When I draw the polygons manualy it works.
Here you can see the differnce:

The first one is my uploaded shapefile and the second one is the manualy drawn shape.

My end goal: To change my shapefile into multipolygon.

*I dont add any code because I don't have code for this and I haven't found any way to do this,yet.


texts: Can you restrict which OSM tags are returned by Overpass API?

texts:
I've been reading up on Overpass syntax and I think I get the basics. For example, let's say I want to get the IDs and coordinates of every town in Andorra that OSM has a record of:

http://overpass.osm.rambler.ru/cgi/interpreter?data=[out:json];area[%22ISO3166-1%22=%22AD%22][admin_level=2];(node[%22place%22=%22town%22](area););out%20ids%20center;

This is great for keeping the data stream light -- say, if you're targeting smartphones and are trying to be conscious of data plan costs.

Now, what if I just want to get the "name" tag of each of these towns in addition to the ID and center? I know how to get all of the tags -- just add " tags" right before the final semicolon -- but that gives me a whole bunch of tags I don't want and a bloaty data stream.

Does anyone know whether/how one can selectively include tags in the Overpass output stream?

That is, not filtering nodes by tag but rather filtering tags by tag name.


texts: Masking a NetCDF file using a shapefile of points with rioxarray

texts:
I took the proposed solution from this question (How to mask NetCDF time series data from a shapefile in Python?) and tried to implement it. However, whenever I apply this solution it seems to set all of my values to NaNs and doesn't mask. Anyone know what's going on here? All of the data details below.
#load in precipitation data
data=xr.open_dataset('/Volumes/Ext HDD 1/Python_data/ERA_precip/data/SA_last10yrs_oct1_nov11/adaptor.mars.internal-1605293227.6880476-7726-17-60b7b292-7985-46a5-ac6e-d3ad5469e87a.nc')

I perform some calculations and up with this percent of normal dataarray that only contains latitude and longitude. It looks like this (pon):

Now, I have a shapefile, which is a list of points in southern Brazil. It looks like this and contains about 2500 points (south_bra_shape):
0       POINT (-50.70833 -26.20833)
1       POINT (-50.79167 -26.20833)
2       POINT (-50.79167 -26.12500)
3       POINT (-54.20833 -25.20833)
4       POINT (-50.95833 -26.20833)

Then I try the solution offered by the previous post:
pon.rio.set_spatial_dims(x_dim="longitude", y_dim="latitude", inplace=True)
pon.rio.write_crs("epsg:32663", inplace=True)
south_bra_shape=gpd.read_file('/Users/eli.turaskyriskpulse.com/Documents/shapefiles/brazil/SC_RGDS_PAR_brazil_soybeans/SC_RGDS_PAR_brazil_soybeans.shp',crs="epsg:32663")
clipped=pon.rio.clip(south_bra_shape.geometry.apply(mapping), south_bra_shape.crs, drop=False)

However, here is the result of clipped:
I do not understand what is happening here or what I am doing wrong.


texts: Get UTM Zone from GeoTIFF (geotools java)

texts:
I've loaded a GeoTiff successfully in with geotools (21-SNAPSHOT), but I'm having trouble accessing the referenced UTM zone.

The code to read in the file is here:

File tiffFile = new File(filePath);
GeoTiffFormat format = new GeoTiffFormat();
Hints hint = new Hints();
hint.put(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
tiffReader = format.getReader(tiffFile, hint);
GridCoverage2D grid = tiffReader.read(null);

I'm inspecting the GridCoverage2D object in my debugger and have found the information in a few places.

The information I'm looking for here is the "UTM zone 16N". Perhaps the zone is given in numerical form somewhere else in the object?

This is inside GridCoverage2D->GridGeometry2D->CoordinateReferenceSystem2D

But for the life of me I cannot figure out how to access this information.
There's no reference to the underlying transformation conversionFromBase in the top level object. None of the accessible methods get me anywhere closer.

Is there an encoded version somewhere in the CoordinateReferenceSystem metadata?

EDIT

The output of the .toWKT() method on my CoordinateReferenceSystem is below

PROJCS["unnamed",
  GEOGCS["GCS Name = unknown",
    DATUM["Datum = unknown",
      SPHEROID["Ellipsoid = unretrievable - using WGS84", 6378137.0, 298.257223563]],
    PRIMEM["Greenwich", 0.0],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH]],
  PROJECTION["Transverse_Mercator", AUTHORITY["EPSG","9807"]],
  PARAMETER["central_meridian", -87.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["scale_factor", 0.9996],
  PARAMETER["false_easting", 500000.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH]]


texts: Split vector by field name and rename output files according to field value

texts:
I have a shapefile with several features and want to split the layer into separate shapefiles. I looked into similar questions recommending Split Vector Layer. However, I could not find a way to specify the name of the output files depending on the value in a certain field.

For example, in the image below, I want to use the field name as the output filenames (like tawitawi.shp for Tawitawi, igat.shp for Igat, etc.), and skip saving the files having blank/NULL values for the name field:

Is there a way to do this, preferably using command line tools to speed up the processing?


texts: Rendering layer in ArcMap using ArcObjects and Python?

texts: I succesfully wrote a script to render a Uniquevaluerenderer in a layer in ArcMAP using symbols from 4 parameters (values) read from the attribute table of the layer. I used ArcObjects and VBA, but with VBA support being discontinued in ArcGIS 10 I now want to convert this script using Python. Can anyone suggest some pointers or sample code.


texts: Creating LAT / LONG grid labels with Northing/Southing and Easting/Westing suffix in QGIS?

texts:
I´m using QGIS 3.6.
I have my map as you can see below:

I would like to have the labels in LAT / LONG coordinates like (40E ; -10S).
Can you please tell me how can I do this?
I saw:
https://www.qgistutorials.com/en/docs/making_a_map.html
Creating lat/lon grid in QGIS?
But it does not solve my issue.


texts: Divide a large polygon into 4 smaller ones having same area

texts: Is there anyway to divide a large polygon into 4 smaller polygons having same area. for example, I have an area of 1 hectare but I would like to divide it into 4 plots having 2500 m2 for each. if anyone knows the solution, please help.


texts: Dissolve adjacent polygons into collections of larger polygons, rather than a single überpolygon

texts:
I have a county's worth of census tracts, and I need to combine adjacent tracts into larger polygons. Here's an example of what I'm looking for: This is four smaller tracts combined into a single large polygon.

I have used Dissolve, but that outputs a single polygon for the entire county which is not helpful for my purpose.

My preferred solution would take a county's worth of census tracts, say 500 or so, and combine them into 200 polygons, without gerrymandering any single group of tracts.

This could be done by hand but I will be doing this for multiple counties.

I'm using ArcMap 10.4.1


texts: Using OpenLayers to work with REST from ArcGIS ImageServer

texts:
I already have ArcGIS Image Server for publishing raster data to my users. I have a third party application that will allow me to add a custom web map as a base layer using OpenLayers2.

I have found endless examples of ArcGIS REST MapServer and WMS services being used with OpenLayers but nothing for Image Server.

Does anyone know ArcGIS REST Image Server support for Openlayers 2 and how to achieve using a REST Image Server with OpenLayers?


texts:
I get the return message NON-VALID PICTURE URL in the POPUP in the CARTO map. when editing in CARTO I get the "Invalid or nonexistent map configuration token". Please forgive my ignorance as I am a newbie at Carto really.
I have followed all the docs and looked at questions but I still can't get the image in the header to display in a popup window.  The image is in the _URL column and is a valid .jpg that I can display by just pasting the URL for it into the address bar of the browser. Is there a size restriction or something like that? https://photos.app.goo.gl/3h6gHXIWQXOW8slo2 is the link to the photo held in Google photos.


texts: Extracting multiple values from list and storing them in data frame using R?

texts:
I have a list of 3 SpatialPolygonsDataFrames. Every element has a number which is the perimeter. I want through a loop to store the numbers into one cell of a data frame like: 154,526,442. I tried to do that like this (in R):

for (i in 1:3){
  if (i ==1){
   datafr$perimeter[1]=mylist[[i]]@data$perim
  }else{
    datafr$perimeter[1]=mylist[[1]]@data$perim;mylist[[i]]@data$perim

  }
}

from this I only get the perimeter of the first element.

How can I store all my values without loosing the previous one every time?(e.g perimeter 1355,1211,..)


texts: ModelBuilder Iteration - Don't want to re-iterate through certain processes?

texts:
I can't seem to figure out how to not re-iterate through the beginning part of my model whereby I import a table from excel into a dbf and join it to the feature classes being iterated. I am copying over data from the imported table by carrying out a join and then just field calc-ing the joined fields over to the appropriate standard fields. The model works but with each iteration it re-imports the excel table, causing the whole process to take longer than needed.

How do I set up the model so that it only imports the excel table once?

Also, is there an easier way to import data from excel and populate into the appropriate fields than the way I'm doing it?

I've attached a screen shot below:


texts: Using Raster Tracer in QGIS v3.6.1 and it's not working

texts:
I've been trying to use raster tracer to trace some lines on simple rasters in QGIS v3.6.1, but when I use the plugin it returns the following error:
2020-06-22T11:48:25     WARNING    Traceback (most recent call last):
File "AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\raster_tracer\pointtool.py", line 217, in canvasReleaseEvent
i1, j1 = self.to_indexes(x1, y1)
AttributeError: 'PointTool' object has no attribute 'to_indexes'
Does anybody know how I can rectify this to make the plugin work or can anyone provide more information about the plugin, as there is very little information on it?


texts:
I'm adding a vector layer that is referencing a mapbox:// url that the user has created via an import. As such, I will not know in advance what kind of features are inside that import e.g. polgyons, lines, points, etc.
How can I apply a smart paint settings that will show the features on the map regardless of type? If you look at studio.mapbox.com they are accomplishing this somehow... just can't figure out how.


texts:
I have a group of about 10,000 POI that are all on streets in the same metropolitan area. I'm about to go through and information about each street from google maps (traffic lane widths).

Is there any way to quickly group POIs into the same street so that I can enter lane widths once for all point on a given street?


texts:
I'm defining "address" as "a physical location, usually visibly numbered".  Numbered was used in the definition to exclude purely geographical features not usually considered to have or need an address.  Visibly was used in the definition to exclude entities that might be numerically coded (like census features) but the number is not typically apparent nor used for directional or delivery purposes.  Usually was used in the definition because exceptions abound (such residences in Carmel).

Enumerating the possible kinds of addresses resulted in:

buildings.
entire floors of buildings.
apartments/offices/rooms in buildings.
cubicles/mailstops in buildings.
loading docks.
slips on boat docks.
slots in trailer parks.
post office boxes and private mail drops.
military APO/FPO addresses.
jail/prison cells.
campground lots.
outbuildings (not usually visibly numbered, but often desireable in delivery databases).
undeveloped parcels (not usually visibly numbered, but often required in address databases).

What other kinds of "physical locations, usually visibly numbered" have I forgotten to include in the above list?  Two questionables are:

hospital beds.
cemetary plots.

Update from comments:

The enumeration (ie, determining how many codes are necessary) is desired to allow coding address kinds to differentiate records in an address database.

Eg, SELECT * FROM ADR WHERE KIND IN (1,3,5).

"Usually visually numbered" (alphanumerics are OK) is preferred to just "coded" as described in the original question. If the location isn't visually numbered, it's unlikely to be the target of directions/deliveries, and would not be kept in the address database.

Just think of the question as "how many different kinds of places do we physically stick a number on to make directions/deliveries easier?"


texts: inf/-inf introduced when opening nc files with raster?

texts:
I am currently working with .nc files and have been struggling - I followed several tutorials and looked into the questions & answers here, but there remains one problem:

I can open and use (/plot etc) the .nc files with raster, however the min and max values in my new raster layer get set to inf and -inf, which leads to problems in my following steps.

When using this:

library(rgdal)
library(ncdf4)
library(raster)

example <- raster("myfile.nc", varname="FAPAR")

I get a "fully functioning" Formal Class Rasterlayer, however when looking at

example@data@min
example@data@max

I get "Inf"/"-Inf"...

Does someone know why and "when" this happens? And what to do?


texts: ArcGIS Geocoding: Error trying to process this table

texts:
I'm using the new what3words geocoding service for ArcGIS Desktop (10.0 SP 5) and keep getting the following error message: "there was an error trying to process this table"

Checked through the threads and nothing applies to my data thus far. it's sample data of a small size, less than 100. The data was provided by what3words, so i know it works! There are no special characters, null values, extra spaces or anything else that typically might cause an issue with geocoding. Tried processing as csv, xls, dbf. Still no luck.

So I know the data is good. I know what3words' geocoding service works (the other features work fine). I know the issue is ArcGIS Desktop reading the data!

A sample of the data (from csv) is provided below.

  what3words Lat Long

  zones.runs.newly   51.549307   -0.440708
  voting.flames.hunter   32.687406   -117.191119
  index.home.raft    51.521251   -0.203586
  shoelaces.upholding.spike  51.054859   0.016843
  loops.keener.cooks 50.866066   2.903226
  spoken.rally.parent    32.622561   -117.070773
  burn.fresh.atomic  32.721688   -116.993317
  timeless.goal.only 51.050655   2.396777
  parent.verge.pyramid   50.617925   3.050013
  operated.comic.abruptly    50.952067   1.913351
  guzzle.digit.switch    50.824668   2.929794
  ducks.august.gymnasium 51.110729   1.259901
  volcano.sorters.tribune    50.865904   2.903397
  cards.punchy.orbit 50.980366   -1.395319


texts:
I run into ST_ColorMap problem in postgis.

Here is the code fragment:

r22=> WITH  sfpoly  AS ( SELECT ST_Transform(ST_GeomFromText('POLYGON(( 3727069.4991851943 4256930.979258035 , 4002242.8010118287 4256930.979258035 , 4002242.8010118287 4440379.847142459 , 3727069.4991851943 4440379.847142459 , 3727069.4991851943 4256930.979258035 ))',3857),4326) as geom), r AS ( SELECT ST_Union(srtm.rast)AS rast FROM o_4_rdepivm1 as srtm JOIN sfpoly sf ON ST_Intersects(srtm.rast::geometry, sf.geom) ),  cx AS (SELECT ST_AsRaster(sf.geom, r.rast) AS rast FROM sfpoly sf CROSS JOIN r )  SELECT  ST_AsPNG(ST_ColorMap(ST_Clip(r.rast, sf.geom), ' 100% 190 185 135  50% 240 250 150  ')) AS rast FROM r CROSS JOIN cx CROSS JOIN sfpoly sf where r.rast is not null;
NOTICE:  More than five elements in colormap entry 1. Using at most five elements
CONTEXT:  PL/pgSQL function st_colormap(raster,integer,text,text) line 85 at RETURN
NOTICE:  Method INTERPOLATE requires at least two non-NODATA colormap entries. Using NEAREST instead
CONTEXT:  PL/pgSQL function st_colormap(raster,integer,text,text) line 85 at RETURN

Can you tell me why?


texts: Setting label colour based on a field value?

texts:
I have created a shapefile with lines which I use to place labels exactly where I want them relative to map features. The shapefile has a field for the name that will appear as the label and also a field called Terrain Type.

I am trying to get the label colour to set based upon the value in the Terrain Type field.

(I am new to QGIS, not using real world data and haven't programmed since using OpenScript with Toolbook about 20 years ago)

Here is what I was trying to do...  Based on another help response I found online, I right clicked on the drop down for the colour in the labels section of the Text section and wrote the following code:

CASE
WHEN "Terrain Type" = "Coniferous  Forest"
THEN Color_rgb (18,95,27)
WHEN "Terrain Type" = "Hill"
THEN Color_rgb (231,238,239)
WHEN "Terrain Type" = "Coastal Plain"
THEN Color_rgb (35,209,16)
WHEN "Terrain Type" = "Bay"
THEN Color_rgb (101,78,240)
WHEN "Terrain Type" = "Scrub"
THEN Color_rgb (155,222,61)
END

It returns an error message saying "Eval Error: Column 'Coniferous Forest' not found".

I have checked the set of similar questions and they don't seem to get me where I need to go.


texts: osm_id in a shapefile vs. the uid and id attributes in an OSM file

texts:
I'm parsing an OSM file using the Python package osmread. The entities parsed have the following attributes id, version, changeset, timestamp, uid, and tags.

On the other hand, I have a shapefile, and I use the Python package pyshp to read it. I then get the following fields from this (road) shapefile:

    [('DeletionFlag', 'C', 1, 0),
    ['osm_id', 'N', 11, 0],
    ['name', 'C', 48, 0],
    ['ref', 'C', 16, 0],
    ['type', 'C', 16, 0],
    ['oneway', 'N', 1, 0],
    ['bridge', 'N', 1, 0],
    ['tunnel', 'N', 1, 0],
    ['maxspeed', 'N', 3, 0]]

My question is, if the road shapefile was obtained from the OSM file, which of the two attributes, id or uid, in the OSM file corresponds to the osm_id field in the shapefile?


texts: Determining the longest side of polygon in QGIS?

texts:
Is there a possibility to determine the longest side of a polygon with QGIS?
As I need to do that for hundreds of polygons I would prefer an automized method (like a plugin).
So far my idea was:

transforming a polygon into a line
cut it on the vertices
measure the length
here I became stuck, how to become the longest line without checking the whole bunch of lines one by one?


texts: Updating SRID (Update_GeometrySRID) does not alter polygons (when it should)

texts:
I am trying to change the SRID of polygons layer, but I'm having trouble visualizing the result.

I the SRID is misreported the shapefile I used to import into PostGIS. The metadata says the reference system is SIRGAS2000(srid 4674). However looking how the data is "drifted" (see image) I suspect the true SRID could be 4618 (SAD69)

So I tried updating the SRID (without changing the underlying points, etc) in Postgis:

CREATE TABLE setor_censitario_2 as select * from setor_censitario

SELECT UpdateGeometrySRID('setor_censitario_2','geom',4618)

and compared the results of the SRID update polygons to to the original polygons in Qgis.

Nothing changed, both layers seem identical.

This puzzled me. Shouldn't they have changed? The parameters of the two "data" (SIRGAS2000 and SAD69) are indeed different in a non-negligible way. So, just to make sure the change of SRID was not supposed to be innocuous, I ST_Tranformed the geom:

create table setor_censitario_4 as
  select geom, ST_Transform(geom,4291) as geom4291Transf, ST_Transform(geom,4618) as geom4618Tranf
    from setor_censitario

now the new geom is different from the original one. ST_Transform, of course, is not the correction I want (in fact I doubled the magnitude of the problem).

So why is Update_GeometrySRID not changing the location of the polygons in Qgis?

The shapefile is available from:
ftp://geoftp.ibge.gov.br/malhas_digitais/censo_2010/setores_censitarios/ba.zip
after unziping I imported with:

shp2pgsql -c -s -I -W LATIN1 file public schema.table ...

UPDATE:

following the [UpdateGeometrySRID help]  (http://postgis.org/docs/UpdateGeometrySRID.html) I also tryed:

CREATE TABLE setor_censitario_3 as select * from setor_censitario
ALTER TABLE setor_censitario_3
  ALTER COLUMN geom TYPE geometry(MULTILINESTRING, 4326)
     USING ST_SetSRID(geom,4326);

results are the same (polygons in Qgis do not move)

UPTADE 2:

I also tried updating the SRID directly in PostGIS, but that did not change the position of the polygons.

To reproduce the situation directly in Ggis, download this file:
ftp://geoftp.ibge.gov.br/malhas_digitais/censo_2010/setores_censitarios/ba.zip
unzip and open 29SEE250GC_SIR.shp in Qgis.

Then add an Open Street Maps Streets layer to the background (I haven't done it but I think it can be done, either by Vector>Open Street Map > Download Data of with the Open Layers Plugin). Or a satellite image from google (also Open Layers).

Zooming to around this location: lon -38.47761, lat -13.01077, Scale 1:8000 you should see something similar to this image:

Green: polygons I want to correct. The perimeters of those should be defined by some of streets (polygons can include more than one city block though).
Red: underlying OSM streets (which I assume are correct)

notice how the green layer is dislocated "northeast" from the red one.

UPDATE 3: changing the SRID to other projections does alter the polygons position. For example, changing the SRIDs from 4674 (Sirgas 2000) to 4225 (Córrego Alegre, also an older Brazillian projection) shifts the polygons south some 50m. It is only in the particular case of 4674 to 4618 that they don´t move.


texts:
I must admit that I'm new when it comes to database administration, but I'm a little confused at the process of my getting my .osm.pbf file into tilemill.

I've run the osm2pgsql process on a .osm.pbf file, it ran, and said that it was completed.  I'm pretty sure that this file goes into PostGreSQL, but I'm unsure of how to access this.

Maybe my understanding of the process is completely skewed?

Not sure, but at this point I know that when I use "psql" in the terminal I'm able to use SQL commands, although I haven't set up any specific account information.


texts: Open GPKG, query layers and save to new GPKG

texts:
I have recently started to use SQL and geopackages queries in R. I know how to open it, access the individual layers (geometry and tables), and select and process individual layers. But, I don't know how to finally store the all the data data again within the new geopackage?

This is the simple workflow:

I'm very new in SQL, so I control my data storage via ArcGIS catalog. And I see, that I can succesfully attach a geometry file, but not the associated table.

Please, what is a correct way to save my geometry, as well as selected table to a new geopackage?

Dummy data are available here: https://www.dropbox.com/sh/zl85ym8wzvzwskd/AAAhODNM5O1lRLb7hOsYCg6za?dl=0

My code:

setwd("C:/Users/.../Dropbox/data")

library(rgdal)
library(dplyr)  # for SQL: http://jsta.github.io/2016/07/14/geopackage-r.html
library(RSQLite)

# Define inputs
gpkgName = "MV_K3243G.gpkg"

# SElect tables to store in new .gpkg
stratumName = "treestratum"
standName = "stand"

# Define new .gpkg
outGpkgName = "out_MV_K3243G.gpkg"

# Read feature geometry and attribute table
fc.shape = readOGR(gpkgName, layer = standName)

# Create connection to SQL database to read table I want to copy
con <- DBI::dbConnect(RSQLite::SQLite(),
                      dbname = gpkgName)

dbListTables(con)  # List and plot all tables in SQL database

# Perform some data.frame operation on the table
stratum_df<-as.data.frame(tbl(con, stratumName))

# save geometry into new .gpkg DOES NOT WORK!!!
writeOGR(fc.shape,
         dsn = outGpkgName,
         layer = "stand",
         driver = "GPKG",
         overwrite_layer=TRUE)

# Write the table to .gpkg
outCon = dbConnect(SQLite(),
                   dbname=paste0(".", "/", outGpkgName))

dbWriteTable(conn = outCon,       # This one does not "stratum_df" table to my .gpkg, only keep the feature geometry "stand"
             name = stratumName,  # How to rewrite it?
             value = stratum_df)

Between data input and data output I would like clip some features and query my data frames, so I don't just simply wont to create a copy of my input geopackage.


texts: Combining two rasters into heatmap using QGIS?

texts: How to combine two rasters so that the heatmap is following according to the relief or shade layer (so that it does not look patched)?


texts: I am unable to run gs:import in WPS Request Builder

texts: I am using the GeoServer version but in that I want to use subprocess for that I have used gs:import but the subprocess is not working. It is giving error while in the other version of GeoServer 2.8. It is working perfectly. I can not get this issue that why my GeoServer is not able to perform the subprocess.


texts: Overlapping transparent polygons won't darken

texts:
A client of ours wanted a map showing the 5-mile and 10-mile radius around business locations today and specifically requested that the overlapping transparent buffer radii grow darker with each other buffer they overlap.  However, after using the buffer wizard and a buffer creator with the business analyst module, the overlapping buffers simply merge together into one shape (similar to performing a Dissolve on them) or, if an outline for them is created, they simply appear to "stack" on top of each other.  This problem did not always occur, does anyone know of a way to cause overlapping transparent polygons to grow more opaque as they overlap?

The software is ArcGIS 10.2.  I believe the problem, after really thinking about it some more, is getting overlapping buffers from point data created by, for example, Buffer Wizard (which creates several buffers as one single layer), to darken as they overlap each other on a map.  If you create each layer individually, then they overlap and darken just fine, but this method is not as reliable.


texts: Problem getting map using XML file with gdal_translate

texts:
The gdal_translate program has an option to define an XML file with the necessary information to download an image using various web map server specifications.  For some reason, the gdal_translate/XML option is not downloading the requested image but going to the web mapping site directly works.  For my needs, I need to go through gdal_translate and use the XML file instead of getting the image directly from the web site.

For example, I have the following XML file that I'm using to download an image:

<GDAL_WMS>
     <Service name="AGS">
         <ServerUrl>https://raster.nationalmap.gov/arcgis/rest/services/LandCover/USGS_EROS_LandCover_NLCD/MapServer</ServerUrl>
         <ImageFormat>tiff</ImageFormat>
         <Layers>show:6</Layers>
         <BBoxOrder>xyXY</BBoxOrder>
        <SRS>4326</SRS>
     </Service>
     <DataWindow>
         <UpperLeftX>-112.2624206542975</UpperLeftX>
         <UpperLeftY>40.5898131802071</UpperLeftY>
         <LowerRightX>-111.0745239257810</LowerRightX>
         <LowerRightY>40.0400730420049</LowerRightY>
         <SizeX>1609</SizeX>
         <SizeY>745</SizeY>
     </DataWindow>
</GDAL_WMS>

I save this file as "test5.xml".  Then I try to run the following command:

gdal_translate -of GTiff "F:\temp\test5.xml" "f:\temp\p37.nlcd.tif"

I get the following output when running this command and the requested image is not downloaded:

Input file size is 1609, 745
0ERROR 1: GDALWMS: Unable to download block 0, 0.
  URL: https://raster.nationalmap.gov/arcgis/rest/services/LandCover/USGS_EROS_LandCover_NLCD/MapServer/export?f=image&bbox=-112.26242065,40.04007304,-111.50641926,40.58981318&size=1024,745&dpi=&imageSR=4326&bboxSR=4326&format=tiff&layerdefs=&layers=show:6&transparent=false&time=&layerTimeOptions=&dynamicLayers=
  HTTP status code: 0, error: SSL certificate problem: self signed certificate in certificate chain.
Add the HTTP status code to <ZeroBlockHttpCodes> to ignore that error (see http://www.gdal.org/frmt_wms.html).
ERROR 1: F:\temp\test5.xml, band 1: IReadBlock failed at X offset 0, Y offset 0
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 0

However, if I go to the following web site, I am able to see the image (this same URL is shown in the output of gdal_translate):

https://raster.nationalmap.gov/arcgis/rest/services/LandCover/USGS_EROS_LandCover_NLCD/MapServer/export?f=image&bbox=-112.26242065,40.04007304,-111.50641926,40.58981318&size=1024,745&dpi=&imageSR=4326&bboxSR=4326&format=tiff&layerdefs=&layers=show:6&transparent=false&time=&layerTimeOptions=&dynamicLayers=

Any ideas about how to fix the XML or work around this problem?


texts: Updating layer extent after changing source in PyQGIS?

texts:
I am running a python script in QGIS 3.8.3-Zanzibar where I am changing the source of a GeoJson file.

I am calling layer.reload() but subsequent calls to layer.extent() still return the old extent.  I can open the attribute table and zoom to the location and see the selected outline.

But the actual layer symbol does not display.

Things that I've tried:

I turned off canvas caching to rule that out.
I have tried these statements (in a variety of combinations):

layer.reload()
 layer.triggerRepaint()
 layer.setDataSource( layer.source(), layer.name(), layer.providerType() )

I figure I'm missing something easy.


texts: converting/parsing geojson into com.vividsolutions.jts geometry

texts:
Does any has any experience with converting string geojson into com.vividsolutions.jts geometry object (point, polygon, etc). Maybe there is some java library I cannot find.

Greets!


texts: Establishing Default Query layer for Google Earth in "Map to KML" export from ArcMap?

texts:
I export a lot of maps for customers who request them from our office. Recently I have seen more and more requests for us to export our GIS into the .KML or .KMZ file type so they can be viewed on Google Earth. I have no trouble exporting my maps into this format using ArcMap's 'Map to KML' conversion tool in the Conversion Tools toolbox.

The issue I'm having is that I get a lot of feedback from customers that they can't properly query the main layer of a map in Google Earth. Instead, Google Earth seems to choose the layer with the least amount of data as the query layer.

I've tried this myself and I can see what they are encountering. I know that turning off the conflicting layers will eventually force Google Earth to query the only active layer, but is there anything that I can do on my end to save them the trouble? This would also in turn save me time from having to troubleshoot this with them when they call. I've attached some screenshots...

I recently created this map using ArcMap and exported it into a .KMZ using the 'Map to KML' tool. The public that we serve loves this as they are able to utilize the strength of GIS in a free program like Google Earth, however when they go to query the features, they get this...

If they don't know much about GIS they can't tell that it's because Google Earth has queried an almost empty layer instead of our main parcels layer. If I turn off the layer Google has chosen to query, it sometimes toggles to some other layer and sometimes goes to the right one.

So far this seems random, but I need to know if anyone knows how to set a priority to the layers when I create the KML on my end so that the finished product is more user friendly for the end users.