texts: How to set up a new gdbase and location in rgrass7 in R, through link2GI, different from the temporary one

texts:
I am using GRASS on Win 10 through the R environment (RStudio, version 1.0.143) with the link2GI package.
I would like to assign a new gisdbase and a new location to my project.

I defined my raster(which has a projection and a datum). Then I called up through
link2GI my GRASS version (7.2):

library(link2GI)
Lammert_Duens_cropped<-raster(paste0(path_output, "Lammert_cropped.tif"))
link2GI::linkGRASS7(Lammert_Duens_cropped, c("C:\\OSGeo4W64", "grass-7.2.1", "osgeo4W"))

It shows, that it has assigned a temporary gisbase and location:

GDAL 2.2.0, released 2017/04/28
gisdbase    C:/Users/KELTOS~1/AppData/Local/Temp/RtmpcbYMjq
location    file1a6812db316e
mapset      PERMANENT
rows        682
columns     423
north       5611931
south       5611249
west        471482.5
east        471905.5
nsres       1
ewres       1
projection  +proj=utm +no_defs +zone=32 +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000 +to_meter=1

So far ok. But if I want to say:

linkGRASS7(Lammert_Duens_cropped,
       gisdbase = "D:/DAS_PROJEKT/Projektarbeit_R_code/R_chaeology/GRASS",
       location = "Init_project",
       gisdbase_exist = FALSE)

R tells me:

  Error in linkGRASS7(Lammert_Duens_cropped, gisdbase = "D:/DAS_PROJEKT/Projektarbeit_R_code/R_chaeology/GRASS",  :    unused
  arguments (gisdbase =
  "D:/DAS_PROJEKT/Projektarbeit_R_code/R_chaeology/GRASS", location =
  "Init_project", gisdbase_exist = FALSE)

And if I say:

link2GI::linkGRASS7(Lammert_Duens_cropped,
                default_GRASS7 = c("C:\\OSGeo4W64", "grass-7.2.1", "osgeo4W"),
                gisdbase = "D:/DAS_PROJEKT/Projektarbeit_R_code/R_chaeology/GRASS",
                location = "Init_project",
                gisdbase_exist = FALSE)

R says the same. R finds Grass, but it doesn't want to assign a new gibsdbase and location, even if I say that gisbas_exist = FALSE.

I know there was an answer (https://gis.stackexchange.com/a/232674/98929) that one can change the default dbase and location, but what am I missing here?


texts: Simplifying variable (if condition) with field calculator in QGIS

texts:
I want to simplify the variable ("SIOSE_CODE") which has values like this:

A(90MTRfr_05PST_05SDN)
A(65PST_35MTRfr)
A(70FDCfr_10PST_10MTRfr_05MTR_05SDN)

The first number (here 90 resp. 65, resp. 70) is a percentage of the land-use class which is most frequently at this location.

Now I want to extract the land-use class (MTRfr, PST or FDCfr), but only if the percentage is higher than 60%.

I tried it with this code:

First I created a new variable "Percentage" and filled it with the output of:

CASE
WHEN (substr("SIOSE_CODE", 1, 1) = 'A' )
OR (substr("SIOSE_CODE", 1, 1) = 'I')
OR (substr("SIOSE_CODE", 1, 1) = 'R')
THEN (substr("SIOSE_CODE", 3, 2))
END

Herethen I want to extract the land-use class:

CASE
WHEN ("Percentage" > 55)
THEN (substr("SIOSE_CODE", 5, 100))
ELSE (regexp_replace("SIOSE_CODE", '[(](.*)?[)]', ''))
END

But the output of my new variable is NULL.

The output of the variable "Percentage" seems to be okay.

Does somebody find the mistake?


texts: Difference between ol.proj.getPointResolution() and ol.sphere.getDistance()

texts:
Is there anyone who can explain the difference between ol.proj.getPointResolution() and ol.sphere.getDistance()?

I have tested it, and ol.sphere.getDistance() is a little bit shorter than ol.proj.getPointResolution().
In addition, I corrected the lengths by elevations with the Pythagoras.
Below an GPX example file.

Why is with "sphere" shortener? "sphere" calculate with a radius and not the direct distance, logical there must be greater. ???
The same difference exists everywhere on the map, also much further north, eg in Europe.

Only getLength(): 1322.334212854104

getLength() with pythagorean(): 1322.4583239930157

ol.proj.getPointResolution(): 1322.3120192140298

ol.proj.getPointResolution() with pythagorean(): 1322.436132435775

ol.sphere.getDistance(): 1320.834204523846

ol.sphere.getDistance() with pythagorean(): 1320.9584565960345

GPX Example (near equador) from brouter.de

<?xml version="1.0" encoding="UTF-8"?>
<!-- track-length = 1323 filtered ascend = 0 plain-ascend = -1 cost=4036 energy=.0kwh time=9.9m -->
<gpx
 xmlns="http://www.topografix.com/GPX/1/1"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
 xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd"
 creator="BRouter-1.4.11" version="1.1">
 <trk>
  <name>brouter_trekking_0</name>
  <trkseg>
   <trkpt lon="11.487967" lat="0.332443"><ele>319.0</ele></trkpt>
   <trkpt lon="11.488969" lat="0.332369"><ele>318.0</ele></trkpt>
   <trkpt lon="11.494446" lat="0.331905"><ele>309.0</ele></trkpt>
   <trkpt lon="11.499806" lat="0.331473"><ele>317.0</ele></trkpt>
  </trkseg>
 </trk>
</gpx>

Example what i mean:

1) is simple length() between two coords. (longest).
2) the earth is round, but a direct distance is "through" the earth. (shortest).
3) the earth is round and the distance is on the earths surface. (best).

TEST:

Result (single-calculation sum, in meter):

Near Equator short distance (4 points), manually set elevation (350,100,0,300):

brouter (inside GPX file) = 1323

ol.proj.getPointResolution() = 1322.3120192140298
ol.proj.getPointResolution() + pythagorean() = 1563.4324417597672

ol.sphere.getDistance() = 1320.834204523846
ol.sphere.getDistance() + pythagorean() = 1562.1085409935195

distVincenty() = 1322.2532472090315
distVincenty() + pythagorean() = 1563.379164847775

Near Equator long distance, 440 meter elevation range (min,max):

brouter (inside GPX file) = 385356

ol.proj.getPointResolution() = 384026.1115163949
ol.proj.getPointResolution() + pythagorean() = 384353.99178101047

ol.sphere.getDistance() = 383596.9242847851
ol.sphere.getDistance() + pythagorean() = 383925.16964882094

distVincenty() = 383079.99691655475
distVincenty() + pythagorean() = 383408.65056785266

Europe near Pyrenees long distance, 430 meter elevation range (min,max):

brouter (inside GPX file) = 144047

ol.proj.getPointResolution() = 142398.39352352163
ol.proj.getPointResolution() + pythagorean() = 142514.5599121801

ol.sphere.getDistance() = 142239.2492071166
ol.sphere.getDistance() + pythagorean() = 142355.53984053037

distVincenty() = 142463.4265797505
distVincenty() + pythagorean() = 142579.5329685958

Fazit: The joke is, the first was i notice, that the ol.proj.getPointResolution() is shorter than brouter. Than i test ol.sphere.getDistance(), and this was more shorter. Last i test distVincenty(), and this was once shorter and some longer than the other (Equator shorter, Europe longer), but shorter as brouter. Logical would be for me, though the ol.proj.getPointResolution() was the shortest and the distVincenty() was the longest. Does anyone have an explanation for that?


texts: How to show labels of overlapping points in the same layer in QGIS?

texts:
If there are 2 points on top of each other in QGIS, only one of them will be labeled

How to make it so that both (or more) points show a label?


texts: Missing roads from OpenStreetMap data with osmosis

texts:
I am currently trying to extract road networks from OpenStreetMap dumps downloaded on Geofabrik. I am using the following osmosis (version 0.45) command to extract roads:

osmosis \
--rbf input.osm.pbf \
--tf reject-relations \
--tf accept-ways highway=motorway,motorway_link,trunk,trunk_link,primary,primary_link,secondary,secondary_link,tertiary,tertiary_link,unclassified,residential \
--tf reject-ways highway=construction \
--lp --wb output_roads.osm.pbf

This should return a file containing all car roads except service roads.

However, some (and sometimes a lot of) roads are missing. I have tried this command on 3 input files:

"Midi-Pyrénées" French region: the road network seems almost complete. The attribute table contains 72010 elements.
all of France: a lot of roads are missing. The attribute table contains only 65803 elements (less than a single region). All types of highways are impacted: primary to residential. There is no any obvious pattern for the missing roads.
all of Europe: almost all roads are missing. The attribute table contains 62092 elements (less than all of France...).

I guess there is some limitation I am not aware of, either on the osmosis command or on the output file. What am I doing wrong here?

This may be unrelated but I have tried using the --used-node option in order to reduce the size of the output files, but this always generates an error (org.openstreetmap.osmosis.core.OsmosisRuntimeException: Unable to create object stream writing to temporary file null.). I have no idea why.

Update 29/08/2018: I also tried using osmfilter (input.o5m --keep="highway=" > output_roads.osm), with the same results. I tried osmosis and osmfilter both on Windows 7 and Debian 9.5 (both 64 bits).


texts: ogr: how to use epsg database for RD New projection?

texts:
I'm trying to assign this EPSG projection to my newly created shapefiles using ogr.
The OGR Projections Tutorial says:

  ...any geographic coordinate system in the EPSG database can be set by it's GCS code number if the EPSG database is available

But this doesn't work for me:

spatialReference = ogr.osr.SpatialReference()
spatialReference.SetWellKnownGeogCS("Amersfoort / RD New")

returns 6 at the prompt.

Doesn't work either (also returns 6 at the prompt by the way):

spatialReference.ImportFromEPSG(28992)

This does work (returns 0 at the prompt):

spatialReference.SetWellKnownGeogCS("WGS84")

If I define the proj4-string manually it works too (also returns 0).

proj4_str = "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889..."
spatialReference.ImportFromProj4(proj4_str)

Conclusion: Not all codes from the EPSG database are available. The spatialReference.XXX() -methods return a code for setting the projection successfully or not(?). Has somebody more experience concerning these techniques and can explain what is going on?


texts: Datatype changed to qlonglong from integer when export to shp and SQL

texts:
I have a view that basis on the table 'mytable' with these columns:

The view created as

CREATE or REPLACE VIEW viewbox.myview AS
SELECT mytable.h AS height,
mytable.geom
FROM mytable;

with columns

When trying to export to shp -> sql

shp2pgsql -s 3857 -I -g geometry /data/relief_poly.shp relief_poly > /data/relief_poly.sql

I get

and in SQL column "height" is float8!

I am not really plugged in to the case why my column datatype changed.


texts: Create union of buffers in PostGIS visible in QGIS

texts:
I would like to merge overlapping buffers created in PostGIS and have as many polygons as there are not overlapping merged buffers, that can be viewed in QGIS.
The answer of this thread seem to be the right path, but I do not understand what g is in the query:
CREATE VIEW buffer40units AS

SELECT
   g.path[1] as gid,
   g.geom::geometry(Polygon, 31492) as geom

FROM
   (SELECT
     (ST_Dump(ST_UNION(ST_Buffer(geom, 40)))).*
   FROM point
) as g;


texts:
I'm working in QGIS with two shapefiles: point and line. Points mark nodes of the lines, and each of them have their own attributes. Lines are related to an attribute table too. The goal is, once the lines are exploded between nodes, to join ONE point attributes to each new line section.
How can I join point attributes to each line section?
(better not using any plugins, so I still cannot download them)


texts: Cutting Fishnet OR Create Parcels in ArcGIS for Desktop?

texts:
What I need to do is the following task in ArcGIS (10.2): I downloaded shapefiles of different countries. For one of these countries, e.g. Austria, I need to define, eg. 400 squares, which means that afterwards Austria is separated in 400 (or any other number) of squares.

As far as I can see I can solve that problem with a fishnet. However, I have some problems with that:

I get to define a Fishnet in any number I want to but:

This number also includes an amount of squares which do not overlap Austria,
I don't know how to cut those squares that do not overlap the shapefile of Austria.


texts: I would like to delineate the polygon of a basin from an already delineated stream work. In the example, I want to extract an approximate polygon from the selected stream network in yellow without using a DEM. Any solution in ArcGIS?


texts: Extracting all cities from a state using PostGIS

texts:
I have written a PostGIS query on a database with OSM data of the state of Baden-Württemberg.

SELECT *
FROM public.planet_osm_polygon
WHERE boundary = 'administrative'
  AND admin_level = '8';

So it turns out that not all cities/towns have an admin_level tag of 8. Is there any way we can further constrain the query to filter out only the towns/cities?


texts: Calculate surface elevation from DEM of pit and known volume with which it is being filled

texts:
I am using ArcGIS 10.1.

I have a DEM of a pit.

I have a known volume of material that I want to put in the pit.

I want to know what the elevation of the flat surface of this material will be once it is all in the pit.

There are a lot of tools for volume calculation that ArcGIS offers, but none seem to relate to this.

What springs to mind is to iterate through a set of constant raster values and see which most closely matches my known volume when I use each constant raster with the DEM in the "cut/fill" tool.

Any other ideas appreciated.

Wil


texts: Requesting map tiles from a wms server in epsg:900913 fails

texts:
I'm trying to figure out how to query a WMS server by using latitude/longitude coordinates. The reason is I'm trying to create an iPhone application and overlay custom tiles with data from the WMS server.

In the iPhone world I can only operate on latitude and longitude and as far as I know, this is the equivalent of EPSG:900913 / EPSG:3857 system.

Now the WMS server examples only lists coordinates in EPSG:32633.

Example url:

http://openwms.statkart.no/skwms1/wms.toporaster2?VERSION=1.1.1&SERVICE=WMS&REQUEST=GetMap&SRS=EPSG:32633&BBOX=196436.79,6685163.95,203635.21,6690907.58&LAYERS=topografiskraster&STYLES=&WIDTH=762&HEIGHT=608&FORMAT=image/png&BGCOLOR=0xFFFFFF&TRANSPARENT=FALSE

After hours of searching on Google, I realize I must be missing something!

Can I query the WMS by using latitude/longitude?
How can I convert between EPSG:32633 and the iPhone format (latitude/longitude) EPSG: 3587?

Any help is appreciated!


texts: Cannot get rid of lock on file geodatabase and feature class created in Python script

texts:
I have a python script, (for ArcGIS 10), that creates a new file geodatabase, creates a new feature class in that file geodatabase and then reads data from an SDE feature class to populate the new feature class.  The script runs fine, except that a lock is placed on the file geodatabase that persists after the script is complete.  The database cannot be deleted with ArcCatalog.  Once the new feature class is deleted, then the file geodatabase can be deleted.

In the script I delete the references to all cursors and rows and yet still the lock remains.  Below is the relevant code.  What am I doing wrong?

import arcpy

OutputFolder = arcpy.GetParameterAsText(0)
FGDName = arcpy.GetParameterAsText(1)
StringInputWorkspace = arcpy.GetParameterAsText(2)

SourceFCAddressPoints = StringInputWorkspace + '\\sde.ERLANDSEN.AddressPoints'

## Create File GDB...
Output_gdb = OutputFolder + '\\' + FGDName + '.gdb'
arcpy.CreateFileGDB_management(OutputFolder, FGDName)

## Set current workspace to to file geodatabase
arcpy.env.workspace = Output_gdb
NewAddressPoints = Output_gdb + '\\AddressPoints'

## Create Feature Class...
arcpy.CreateFeatureclass_management(Output_gdb, "AddressPoints", "POINT", "", "DISABLED", "DISABLED", "PROJCS['NAD_1983_StatePlane_Washington_North_FIPS_4601_Feet',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Lambert_Conformal_Conic'],PARAMETER['False_Easting',1640416.666666667],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-120.8333333333333],PARAMETER['Standard_Parallel_1',47.5],PARAMETER['Standard_Parallel_2',48.73333333333333],PARAMETER['Latitude_Of_Origin',47.0],UNIT['Foot_US',0.3048006096012192]];IsHighPrecision", "", "0", "0", "0")

## Add fields
arcpy.AddField_management(NewAddressPoints, "AddressPointID", "LONG", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

curSourceAddresses = arcpy.SearchCursor(SourceFCAddressPoints)
curNewAddressPoints = arcpy.InsertCursor(NewAddressPoints)

for rowSourceAddress in curSourceAddresses:
    rowNewAddress = curNewAddressPoints.newRow()
    rowNewAddress.SHAPE = rowSourceAddress.SHAPE
    rowNewAddress.AddressPointID = rowSourceAddress.OBJECTID

del rowNewAddress
del rowSourceAddress
del curSourceAddresses
del curNewAddressPoints


texts:
I have edited a field into the attribute form of a QGIS project as a "Value Relation" to create a checkbox where I can choose few possibilities. I have a tree sickness and the field I am talking about is the LOCATION where the sickness is visible: Roots, leaves, trunk, branches, twigs. The thing is that the output of all of that when I open my attribute table after typing the data it is like this in case the sickness is in all parts of the tree:
{"Roots","Trunk","Braches","Leaves","Twigs"}

Is there anyway where I can transform this into a normal string when it comes to be readable on the attribute table like Roots, Trunk, Branches, Leaves and Twigs (without the and also good solution but for sure without curly braces and quotation marks)?


texts:
I am starting to learn python.
I have a set of shapefiles (points) and a set of GeoTiff raster files (clipped, and with landcover classes).
Now I want to intersect them. As result I would like to get specific information about each single point. (which land cover class intersects with the point when setting certain buffers e.g: 50 m, 100 m, 150 m)
How to handle such a task using python?


texts: Documenting field names that Network Dataset Wizard looks for to create network attributes?

texts:
Is there a definitive list of the field names that the Network Dataset Wizard looks for when it runs, and creates network attributes from?

For example, the Network Dataset Properties below resulted from taking all the defaults on the way through and seems to have automatically created three network attributes:

Length - from the Shape field
Minutes - because it found a Minutes field
Oneway - because it found a One_Way field

I would have expected these "special" field names to be fully documented at, or on a page linked from, the Creating a network dataset page, but they appear to only be done to this level:

  Be sure your sources have fields that represent your network impedance values—distance, travel time, and so on. For best results,
  name these fields as the units of your impedances, since these fields
  will automatically be detected by the New Network Dataset wizard—for
  example, you might want to name your travel time field Minutes. For
  edge sources, if the impedance values differ based on direction of
  travel, provide a separate field for each direction of travel, for
  example, FT_Minutes and TF_Minutes.
  If you are modeling one-way streets, be sure that your edge sources have a field specifying one-way street information. The New Network
  Dataset wizard recognizes a string field named One_Way or Oneway,
  creating evaluators that interpret its values as follows:

  FT or F indicates a one-way street permitting travel only in the digitized direction of the edge.
  TF or T indicates a one-way street permitting travel only against the digitized direction of the edge.
  N indicates a street that does not permit travel in either direction.
  Any other value indicates a street that permits travel in either direction.


texts: Is there any web, scripts, or tools to calculate number of registered OSM user in each country?

texts:
I will have talk about OSM in next few week and I just wondering if there any tools, scripts, or websites to showing the statistic of number registered OSM User per country.

Because I have an idea to displaying the number of registered user (for example by countries or one country per years) in my talk in order to compare which country that have the biggest number of registered OSM users.

One website that i found is here: https://osmstats.neis-one.org/

But it only showing the whole world, not in specific country. When i moved into the countries tab in that web, it only show the number of active user each day.

Any advice?


texts: How do you change a GraphicsLayer's geometry in the ArcGIS API for Javascript 4.0

texts:
The setup:

Creating a Point with a longitude and latitude
Creating a SimpleMarkerSymbol
Creating a Graphic using the Point and SimpleMarkerSymbol
Adding the Graphic to the view's graphics

What's I'm trying to do:

Change the point's latitude and longitude
Have this change be reflected in the map.

What I've tried:

Just changing the point's latitude and longitude.
Removing the Graphic with the point, then changing the latitude and longitude, then adding the graphic back to the view.
Looking through the ArcGIS API reference for Javascript 4.0 for anything like refresh/update/redraw/... on Map objects, MapView objects, and GraphicsLayer objects

Some sample code to play with:

I took the following from an ArcGIS example and modified it to show as simply as possible what I'm trying to do.

<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no">
  <title>Get started with graphics - 4.0</title>

  <link rel="stylesheet" href="https://js.arcgis.com/4.0/esri/css/main.css">

  <button id="movePoint">Move Point to 0,0</button>


texts: Determine lat/lon intersection on earth from point in space along 3d line

texts:
I have a known point in space with a lat,lon,alt.
It has a camera with a line of sight and a angle of field of view. I can take the central point and calculate each corner vector of the field of view, but..

How do I calculate the intersection with the earth (not necessarily at MSL) from these 3D lines from a point in space? I already know the height of this intersection, but need the lat/lon.

If the earth was flat and square - then that is simple trigonometry.
But curvature of the earth and rounded longitude lines has me stumped; and I can't find anything on the net that seems to provide a first step..
I will implement the solution in c#.

Any advice from the 3D gurus here?!


texts: ArcPy Get Value of Field for Selected Record

texts:
I am building a script to populate a field in one feature class based on a field in another feature class.  First I have a leases layer which is all small polygons.  Then I have an index grid layer which is large polygons.  I want to know which leases exist within which index grid polygon.  The index grid layer has an "Id" field, and I want the Id value from the index grid that each lease exists within to be written to a new column in the leases layer.

For example, index grid polygon with an 'Id' of '7' contains leases '14' '6' and '12'.  The final result should be a column where lease 14, 6, and 12 all have a value of '7' in the new column.

Here is what I have so far:

import arcpy, sys

grid = arcpy.GetParameterAsText(0)
leases = arcpy.GetParameterAsText(1)
workspace = arcpy.GetParameterAsText(2)
new_field = arcpy.GetParameterAsText(3)

if arcpy.ListFields(leases, new_field):
    sys.exit(arcpy.AddError("Error:Field Already exists"))

arcpy.AddField_management(leases, new_field, 'SHORT', '', '', '')

grid_lyr = arcpy.MakeFeatureLayer_management(grid, 'grid_lyr')
leases_lyr = arcpy.MakeFeatureLayer_management(leases, 'leases_lyr')

fields = ['file',new_field]
with arcpy.da.UpdateCursor(leases_lyr,fields) as cur:
    for row in cur:

        leases = row[0]
        where = '"file" = \'{}\''.format(leases)
        arcpy.SelectLayerByAttribute_management(leases_lyr, 'NEW_SELECTION', where)

        arcpy.SelectLayerByLocation_management(grid_lyr, 'CONTAINS', leases_lyr)
        map_no =  #I think this is all that is missing

        row[1] = map_no
        cur.updateRow(row)

The value 'map_no' represents the index grid Id which the lease belongs to.

As I have a comment in my code, I think there is only one line that is missing...in the original script which this is derived from, this line uses the arcpy.GetCount_management, but obviously it needs to be something different.  It just needs to be some kind of get value function, but I simply cannot figure it out.

I derived this work so far from a script that counts the number of points in in a polygon, and then writes that value to a new column in the polygon layer.  I mention this because maybe there is a better solution.

Using ArcMap 10.1.


texts: Using gdalbuildvrt in R gives Warning Message ... had status 1?

texts:
I want to make a mosaic out of several raster files that are next to each other, all including 4 bands. Due to saving space on the server I need to build a vrt and want to write the code in R, so I can process this action more often without too much work. Unfortunately I cannot get rid of the Warning message:

    Warning message:
    running command '"C:\Program Files\QGIS 2.16\bin\gdalbuildvrt.exe"
    -overwrite
    "S:/Raster/Ortho_Mosaic.vrt" "S:/Raster/Input/*.tiff"'
    had status 1

I use R and gdalbuildvrt should use gdal from QGIS 2.16.

My packages are rgdal, raster and gdalUtils.
My code looks like this:

    folder <- 'S:/Raster/'
    gdalbuildvrt(gdalfile = paste(folder, '/Input/*.tiff', sep = ""),
    output.vrt= paste(folder, '/Ortho_Mosaic.vrt', sep = ""),
    intern = TRUE, overwrite = TRUE))

I tried a lot of different ways for my input files:  use an extent, use single files, all .tif files or changed tif to tiff ...
The output file is not produced even though it is only a warning and not an error message. All files are in the same coordinate system; EPSG 31467

Does anyone have an idea what the problem could be?


texts: Masking pixels inside the polygons using Google Earth Engine?

texts:
I am working on a classification problem on agricultural fields. To make things easier, I first identified road, cities, water, fields, and etc. I am planning to eliminate road, cities, and water out of the image to improve my model. I got all these areas as polygons in a .shp file and I want to mask out these areas.

Wherever I check, all applications are for clipping a polygon.

Is there a method to mask the polygons I provided in GEE code editor?


texts:
I have a shapefile with a number of fields.  I have one field (let's call it Field A) that I need to populate based on another field (let's call it Field B).

If Field B is greater then 0 I need it to populate Field A with the text 'Fill'
If Field B is zero then I need it to populate Field A with the text 'NoChange'
If Field B is less then 0 I need it to populate Field A with the text 'Drain'

I should know what the Python code for this needs to be.


texts: Calculating field on insert (real time) using QGIS?

texts:
We use QGIS as a field application. We like to calculate with a input value and then also calculate with the result of the calculation

Example:
Liters: 1000 ->
Liter per container: "this has to be Liters / 3
With the result of Liter per container  we also want to calculate

Did use the field calculator, but this doesnt work real time.

IS there a option to do the calculation real time?


texts: How to divide and split an arc into specific lengths?

texts:
I have an arc shapefile for a city which represents the rail road route. It is a single attribute arc that I would like to separate into segments of 90 ft.

For example, if the standard track length for my railroad is 90 ft and the length of the route is 9000 ft, I would like to have 100 attributes, head to tail, in a single arc but with separate attributes.

Is this possible in QGIS?


texts: When working in QGIS (3.4.6), how do I get a layer symbol(point) to remain its original size? In this case I have a point which enlarges in size when I zoom out, but it should remain it's original size. How does this work?


texts: Uncaught ReferenceError: require is not defined ArcGIS javascript API

texts:
Below is my code.I'm tryng to display a GeoJSON layer from ArcGIS online. I get the error,
Uncaught ReferenceError: require is not defined.

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta http-equiv="X-UA-Compatible" content="IE=7,IE=9">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
    <title>ArcGIS GeoJSON Layer</title>

    <link rel="stylesheet" href="https://js.arcgis.com/3.9/js/esri/css/esri.css">
    <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.1.1/css/bootstrap.min.css" >
    <link rel="stylesheet" href="https://esri.github.io/bootstrap-map-js/src/css/bootstrapmap.css">

Can anyone help me on this?


texts:
Bom dia, estou com alguma dúvidas a respeito do QGIS WEB Cliente que não consigo achar uma solução. Eu sei que tenho que fazer a instalação do Qgis Serve como esta escrito neste tutorial Anita Graser QGIS Server on Windows 7

I also know that I have low ZIP QGIS Web Client, but where should I put it in Windows? After that should follow this documentation

In my attempt put the contents of the ZIP file into C: \ OSGeo4W \ apache \ htdocs \

The map model HelloWorld works

But I have no idea where to go, my project is stopped.


texts:
Looking to create a distance to wetland raster. However this travel distance must be through the water only.

I have tried setting land to a high cost (e.g. 99) however doing so I have to set each water cell to its own "cost" however I don't want to calculate relative cost, I want the actual distance in Metres.


texts: Access to SpatiaLite .dumpshp via Python?

texts:
I'm trying to interact with a SpatiaLite database via Python on my Windows machine.

I've managed to connect to the database and load the SpatiaLite extension, which means that I can pass SQL queries in from my Python code e.g:

from pysqlite2 import dbapi2 as sqlite
myPath = r"C:\SpatiaLite\test-2.3.sqlite"
conn = sqlite.connect(myPath)
conn.enable_load_extension(True)
curs = conn.cursor()
curs.execute("select load_extension('libspatialite-1.dll')")
curs.execute("""SELECT name, peoples, AsText(Geometry) from Towns where
                peoples > 350000 order by peoples DESC""")
# Print the selected records
for row in curs:
    print row

So far so good. I'd now like to know how to export a table as a shapefile. Using the command prompt in spatialite.exe I can do this using .dumpshp, but I can't figure out how to access this command from Python. Any ideas, please?

I was hoping for something like this:

curs.execute("select dumpshp('Towns Geometry shape_towns CP1252 POINT')")

but perhaps it's more complicated than that?

Thanks very much!


texts:
I would like to calculate the difference of area between two polygons. I tried using the Difference tool, but all it does is create lines that are different between my two polygons. I also tried the Symmetrical difference tool, but it for some reason it did not work and just generated the input layer out instead of giving me the extra polygons of my two layers. Are there other ways for me to find the area difference between the two polygons? In my attridutes table, all I have for the two polygons are their area.

And the screen capture of my case.

Here is what I meant. The purple and green polygon are two different layers (Which largely overlap with one another with the green layer being covered by the purple layer). I would like to know the area difference between the two, which is now the non covered parts you can see of the green polygon.

I am new to QGIS so I do not know handy methods that can do.


texts:
I am interested in finding the intersection of a Ray and a Surface in ArcObjects.
I already know of two solutions.  First solution is to use the Locate method of ISurface:
Dim intersectionPoint As IPoint = surface.Locate(ray, hint)

The other solution is to use the GetLineOfSight method of IGeoDatabaseBridge (where the variable point is a calculated point along the path of the ray):
Dim geoDatabaseBridge2 As IGeoDatabaseBridge2 = New GeoDatabaseHelper
geoDatabaseBridge2.GetLineOfSight(surface, baseRay.Origin, point, intersectionPoint, Nothing, Nothing, Nothing, False, False)

Surprisingly, the GetLineOfSight method is much faster than the locate method.  But I am looking for an even faster solution.
So my question is:  Is there any other (and faster) ways to find the intersection between a ray and a surface?  Or can you point me in the direction of literature that explains the calculations I need to roll my own?
Thanks!
UPDATE:
I have started the process of creating my own Line of Sight algorithm.  I have posted my current solution description in a blog post.  Thanks again for all of the suggestions.


texts: Export lines or graphics without losing the details / generate graphics without having lines overlapping

texts:
I am new to ArcGIS.

How do I easily reproduce a spatial diagram such the one attached (Picture); from shapefiles on ArcGIS.

I very often have to reproduce networks on Powerpoint at different scale without loosing details such as number of lines from A to B.
Exporting images directly from ArcGIS at certain scales makes those details difficult to notice, as lines close to each other tend to look like one.

Is there a way to export lines or graphics without losing the details, or is there a way to at least generate graphics without having lines overlapping?

(Note that in reality there only is about 30meters between the parallel lines, and that the distance between the red dot and the green square is 20km - I need to keep the relative location of point and the relative route of the lines, but at the same time represent all the lines although on the ground they are very close to one another )


texts: Calculate distances between polygon centres in QGIS

texts:
Is there a practical way to calculate distances between polygon centroids in QGIS? I have a UK city region map divided into administrative units (1 ward = 1 polygon). I want to represent their population growth in function of the distance to the core city and the best way to do it seems to be looking for the geometrical 'centroid' of each polygon and measure its distance to the 'centroid' of the core city polygon. I can do it manually with the ruler tool but it takes a long time and is not very accurate.
Thanks


texts: QGIS crashes unexpectedly when using Raster Clip tool

texts:
Running QGIS 2.18 on Mac OS 10.14.
Python 2.7 and GDAL 2.1 installed.

Trying to clip raster by polygon (raster--> extraction --> clipper) I get error message: qgis 2.18 clipper The process failed to start. Either the invoked program is missing, or you may have insufficient permissions to invoke the program*.

Then QGIS crashes, gives the following message:

Process:               QGIS [4831]
Path:                  /Applications/QGIS.app/Contents/MacOS/QGIS
Identifier:            org.qgis.qgis2
Version:               2.18.28 (2.18.28 [exported])
Code Type:             X86-64 (Native)
Parent Process:        ??? [1]
Responsible:           QGIS [4831]
User ID:               501
Date/Time:             2019-04-17 23:41:07.988 +0200
OS Version:            Mac OS X 10.14.2 (18C54)
Report Version:        12
Anonymous UUID:        920AA060-0540-B21D-7AD7-A6AF2928F792
Time Awake Since Boot: 4700 seconds
System Integrity Protection: enabled
Crashed Thread:        0  Dispatch queue: com.apple.main-thread
Exception Type:        EXC_BAD_ACCESS (SIGSEGV)
Exception Codes:       KERN_INVALID_ADDRESS at 0x00003e10c0113cd0
Exception Note:        EXC_CORPSE_NOTIFY
Termination Signal:    Segmentation fault: 11
Termination Reason:    Namespace SIGNAL, Code 0xb
Terminating Process:   exc handler [4831]
VM Regions Near 0x3e10c0113cd0:
    JS JIT generated code  0000249fa5400000-0000249fdd400000 [896.0M] rwx/rwx SM=NUL  reserved VM address space (unallocated)
-->
    MALLOC_NANO            0000600000000000-0000600010000000 [256.0M] rw-/rwx SM=PRV

What's going on? Why does raster clipper not work? More lines follow.


texts: OpenLayers wrapDateLine Problems with ArcGIS93Rest

texts:
I'm having issues viewing a ArcGIS93Rest layer on both sides of the dateline -- the data doesn't cross the dateline. The problem is that I only see data for one side of the dateline until I pan west and vice versa (see attached pics at bottom).

I'm trying to rule out that it's an OpenLayers problem. The wrapDateLine option should work with the ArcGIS93Rest layer type because it's a subclass of Layer > HttpRequest > Grid same as OpenLayers WMS Layer. I have no wrapping problems with WMS layers created with UMN MapServer.

However, I can't set the wrapDateLine option as a hash extra option on the ArcGIS93 layer but only after initializing a layer like so. I'm not sure this is right. Is it?:

layer = new OpenLayers.Layer.ArcGIS93Rest('L' + cfg.id, cfg.url + '/export',
    { layers: 'show:' + cfg.layer, transparent: 'true',
        f: 'image', format: 'png', size: '256,256',
        bboxSR: '900913', imageSR: '900913'
    });
layer.setOpacity(cfgOpacity);
layer.wrapDateLine = true;

Questions:

Does anyone know if wrapDateLine works with ArcGISRest93? Have you seen it work?
Has anyone had this problem with ArcGIS Server data before? I'm wondering if it has to do with their service. As I said before, I have no problem with WMS layers from UMN MapServer


texts: How to change the color of a tile layer?

texts:
If I generate a very large number of map tiles (300,000), is there an easy way to let a Google Maps user change their color?

For instance, is it possible to change the color scheme of a map tile layer by overlaying it with other images (red, blue, green) with varying degrees of opacity?

My map tiles are showing household income at the census tract level. So that is a single numeric value represented by a color.  So I don't risk losing information if I add (or remove) red, blue or green from it.


texts: Select by attributes - records with the largest attribute value

texts:
Using 'Select By Attributes' how can I select records from a shapefile database table where a numeric field contains the largest value.
Need to use MAX function or something like that but can't get the syntax right.


texts:
I work on arcview 10.3 and i try to find layer called "ways" in several mxds. i try this code but it does not work. anyone know how to do it?:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.workspace = r"C:\Project"
for mxd in arcpy.ListFiles("*.mxd"):
    print mxd
    mapdoc = arcpy.mapping.MapDocument(r"C:\Project\\" + mxd)
    df = arcpy.mapping.ListDataFrames(mapdoc, "Layers")[0]
    for lyr in arcpy.mapping.ListLayers(mapdoc, "", df):
        if lyr.isGroupLayer == True: continue
        if lyr.symbologyType == "UNIQUE_VALUES":
            if lyr.name == "ways":
                print 'yes'
    mapdoc.save()
del mxd


texts: Clipping shapefiles with different CRS in QGIS 2

texts:
In QGIS 2.18.22 I have to different vector-shapefiles and I would like to clip the overlay.

Unfortunately, the shapefiles did originally not have the same CRS so I saved them in the desired CRS (EPSG:4326, WGS 84). After that, I opened a completely new project (EPSG:4326, WGS 84) and opened the new shapefiles (OTF is turned off).

However, there is still no overlap between both shapefiles, but I can zoom to each of the shapefiles individually. Clipping is of course still not working.

Do you have any ideas on how to solve that problem? I checked the similar questions but I did not find any solution for my problem.


texts: cutting and extending one polygon layer in function to the other layer boundaries

texts: I'd like a methodology that can edit my vector layer, extending and cutting the features in the polygon layer "cap_agricola" within the boundary layer "sp_edr", I would like to solve this problem in ArcGIS, but I can do without restrictions in GRASS or QGIS.


texts:
Im am trying to export my contours from ArcMap to AutoCad from a .tiff. I am able to get the contours extracted and opened in AutoCad but they come out flat. I made sure I was inserting the correct 'contour interval' when extracting contours from the raster but they still come out flat.
Can any one please help out?


texts: How to create an image of maximum pixel values from an input layer-stack image in ERDAS Imagine?

texts: I have 12 MODIS images (one of each month). I want to create a single image having maximum annual pixel values for each pixel. I did layer stacking operation in Erdas Imagine 13, but how to extract the maximum pixel value ?


texts:
I have created a fishnet around a city-plan.

Because the fishnet is huge, the loading time is of course also big, too big to continue working smoothly.

So I'd like to clip the fishnet to only the zone of the city I need:

So I become something like this when zoomed in:

So I assume I'll need the Clip tool for this, but based on what "Clip Feature" can I clip the fishnet?
Probably some kind of square, but how do I sketch my desired square?


texts: Finding points in one layer within range of time of points in another layer using QGIS

texts:
I am very new to QGIS and SQL.
I have two layers, each is populated with points and temporal data. One layer has information on battles occuring in Iraq and Syria, and the other has information on terror attacks. I would like to extract information about the distance of each terror attack to the nearest battle, within 30 days.
I have two SpatiaLite layers and have tried the following query to create a new virtual layer:
SELECT battles.data_id as id_battles, terror.fid as id_terror,
       ROUND(MIN(ST_Distance(acled.geometry, fullset.geometry)),2) AS distance
FROM ACLED AS battles, FullSet_Prio2 AS terror
WHERE battles.'date' BETWEEN terror.'date' AND (terror.'date' - 1 month)
GROUP BY terror.fid
ORDER BY distance DESC

This returns a syntax error at the -1 month line. +
The query works if I use == instead of WHERE...BETWEEN...AND but that's not the information I need.
Could anyone help me figure out the correct syntax for substracting from dates or do I need to do something else instead?
I tried creating a virtual layer using this:
WHERE day(age(fullset.'date', acled.'date')) < 30

This doesn't result in an error, but results in empty rows.


texts: Remove polygons of one layer from another polygon layer in QGIS

texts:
I have two polygon layers, A and B, both in the same projection. I am trying to create a new layer, C, which is all polygons in A less the areas where polygons in A and polygons in B intersect.

I am using v2.14.1 - Essen.

I have checked the geometries are valid in both layers, removed duplicate notes and dealt with crossovers. Running the validity query in QGIS shows that there are no errors.

I ran the 'difference' geoprocessing tool, and, for the most part, this works. C is created, which includes all A polygons, less those areas that overlapped with B polygons.

However, there is one polygon in A which, for some reason, does not seem to work as part of the process. When I create C, the entirety of this polygon has disappeared, despite there only being a small overlap with several polygons in B.

Can anyone point me towards a reason why this might be?

(Apologies for A, B and C, but thought I would generalise the question)


texts: Fetching postgis data through php json and passing to leaflet map

texts:
I am new to this and trying to explore leaflet and geojson a bit. Please help me with this.

This is my php page which I am using to fetch data from PostGIS. (the data opened successfully in the browser)

$sql = "SELECT *, (ST_AsGeoJSON(geom)) AS Geojson FROM
   public.village_boundary_metro";
$rs = pg_query($dbconn4, $sql);
if (!$rs) {
   echo "An SQL error occured.\n";
   exit;
}
$output    = '';
$rowOutput = '';

while ($row = pg_fetch_assoc($rs)) {
  $rowOutput = (strlen($rowOutput) > 0 ? ',' : '') . '{"type": "Feature", "geometry": ' . $row['geojson'] . ', "properties": {';
  $props = '';
  $id    = '';
  foreach ($row as $key => $val) {
      if ($key != "geojson") {
        $props .= (strlen($props) > 0 ? ',' : '') . '"' . $key . '":"' . escapeJsonString($val) . '"';
      }
      if ($key == "id") {
        $id .= ',"id":"' . escapeJsonString($val) . '"';
     }
   }

  $rowOutput .= $props . '}';
  $rowOutput .= $id;
  $rowOutput .= '}';
  $output .= $rowOutput;
}
$output = '{ "type": "FeatureCollection", "features": [ ' . $output . ' ]}';
echo json_encode($output);

And now I try to extract this resultant data through Ajax as a geojson layer on a leaflet powered map.(ajax code follows)

$.ajax({
  type: "GET",
  url: 'geojson_entry_postgis.php',
  data: data,
  dataType: 'json',
    success: function (data) {

  var geojson = new L.geoJson(data).addTo(map);
}})

Also tried this way..

var geojson = L.geojson('geojson_entry_postgis.php?geotable=village_boundary_metro&geomfield=geom',
    { style: {"color":"#ff7800","weight":2}, maxZoom: 6, onEachFeature:
       onEachFeature }).addTo(map);

However the data fetched from postgis does not display as geojson layer on the map as expected. Please let me know where I am going wrong and Am I following the right process??

P.S. Ajax function runs fine in other examples.


texts: Why does destroying a popup cause a "this.events is null" error?

texts:
I am trying to create dynamically generated popups in OpenLayers with a mouseover event. According to the documentation, it is recommended that the popups be destroyed after use to avoid memory leaks. I set up an event listener that destroy the popup upon a mouseout event as follow:

marker.events.register("mouseout", marker, function(event) {
    popup.destroy();
});

While testing on jsFiddle upon map initialisation, I am able to create and destroy the popup on the first mouseover and mouseout. Subsequent attempts of mouseout result in a "this.events is null" error appearing on my javascript console log. It seems to me that the first instance of popup.destroy() might have destroyed more than just the popup itself, causing the second instance to fail.

Can anyone enlighten me on how to resolve this issue without deviating too much from the current setup?


texts: ArcGIS Print Javascript API return syntax error after installing SSL certificate

texts:
I've been facing problems with my printing service using ArcGIS JavaScript API 3.22 and ArcGIS server 10.1 after installing SSL certificate on the ArcGIS server.
The printing service return the SyntaxError: Unexpected token < in the init.js and Print.js file. I've done some digging and found that this problem will occur while printing large amount of graphics on the FeatureLayer. Is there any configuration I need to do on the ArcGIS server after installing SSL certificate so that it can print large amount of graphics because i did not face this problem while printing on ArcGIS server without SSL certificate.


texts: I am trying to merge 4 TIFF files and the resulting merge is coming out with some odd colour distortions. Does anyone know how to resolve this?


texts:
I have an input jpg image (sat thumbnail) and I know a number of gcps for this image. I want to reproject the jpg to have a north up image. This is what I do:
./gdal_translate -of GTiff -a_srs EPSG:4326 -gcp pixel line lon lat -gcp pixel line lon lat ... etc... myinputfile.jpg temp.gtif
./gdalwarp -t_srs EPSG:4326 -r near temp.gtif myoutputfile.tif

The result is "almost" correct, although I do get some strangest behaviours sometime, eg the image folds on itself, etc... very artistic but not good for a GIS application!
If anybody can see anything wrong with the way I proceed or has a better idea, please advise!
Edit
Alright, it is working now. I believe the "weird" results were down to the fact I used a temporary file that was not deleted between gdal calls. As soon as I had the temp file removed (rm temp.gtif) between calls the results were fully normal.


texts:
I'm working with Leaflet and nodejs and when a polygon layer is loaded when the map view is at world map, the polygon don't shows, I have to zoom in at approximately zoom 10 or 13 to show polygons.

I was trying to use this plug-in but realy don't know how to use it. The lib https://github.com/oliverroick/Leaflet.Deflate

To represent the map I use this

    var layers = {
      type: 'FeatureCollection',
      features: features,
      properties: {}
    };

    showLayers = L.geoJson(layers);

    leafletData.getMap().then(function(map) {

      showLayers.addTo(map);
    });

And I try to use it but don't know where?


texts:
Is it:

Authors: John Doe, Jane Doe

or

Author: John Doe, Jane Doe.

It is a small difference, but I wondered how to do it correctly.


texts:
I've created simple model. It worked fine, but next time after opening ArcGIS I saw model in that way:

Used tools are not orange but white and I can't click on them (to edit something).

When I run my model it still works properly, but I'd like to change something which is impossible.

Any ideas what's going on?


texts: How to set feature coordinates to values in attribute table?

texts:
I have some shapefiles that were produced programatically, and everything came out okay EXCEPT that the geometry information was not written correctly. The coordinates for each feature are null. Luckily it also made a duplicate x,y,z column with the coordinate information in there.

How can I rewrite the coordinates for each feature to pull from the attribute data? Particularly in Python/ArcPy so that I could script it (Java/GeoTools can also be an option).

Edit: These are point data (which makes sense as one x-column, one y-column, one z-column).


texts: Buffer not displaying and reference points disappearing in QGIS?

texts:
I have created 'reference points' on QGIS 3.6.3 and ensured I am using a metric coordinate system. The reference points and all data display correctly but when I attempt to add a 150m buffer around these 'reference points' (using either processing toolbox -> vector geometry -> buffer OR processing toolbox -> GDAL->vector geoprocessing -> buffer vectors) the reference points disappear and the buffer does not appear on the map, despite being created as a layer.

Does anyone have a solution to this problem?


texts: Seeking Singapore postal districts boundaries?

texts:
Does anyone know where to get GIS data for Singapore postal boundaries?

I need a shape or KML file which I can then overlay onto google maps to show a neighbourhood break down of Singapore.


texts: Identify if point over land or ocean using lat/lon in R with confidence metric

texts:
How can I determine whether a coordinate is above land or ocean with a confidence metric in R?
I've used the sf package method described in another question: Given Lat and Lon identify if point over land or ocean using R
#Obtain latitude and logitude values.

lat <- filtered.records.coord$lat
lon <- filtered.records.coord$lon

#Place the values into a data frame.

points <- data.frame("Var1" = lon, "Var2" = lat)

#Find which points fall over land

pts <- st_as_sf(points, coords=1:2, crs=4326)
ii <- !is.na(as.numeric(st_intersects(pts, world)))

#Combine GPS results with record information.

filtered.records.coord$land <- ii

However, is there a way to determine if a coordinate is by the coastline for example and not far inland? Or some sort of confidence to the TRUE or FALSE offered by the sf package? I wish to identify borderline points.


texts: How do I avoid text field truncation when saving a table join?

texts:
Using QGIS I've joined tabular text data (all in fields of less than 255 characters) to a shapefile via a CSV file. All looks good in the joined table but when I save to a new shapefile to make the join permanent the long text fields get truncated to between 50-100 characters. How can I avoid this? It seems the problem is in the saving/exporting process since the initial join is fine.

I'm running QGIS 1.7.0 on Mac 10.6.8

Thanks,
Matt


texts: Clustering on Geoserver layer appearing out of the map

texts:
We are using Geoserver and openlayers to display maps. Geosever layer provides list of sites, since there were many we did a clsutering using Styles.

But the problem is when we zoom out completely some of the clustered points are appearing out of the map.

Please help me to know the actual issue with this. i want the clustered point to be within the bounds of the map.

go to aekos.org.au and select user defined locations. Add some layers to see what I am talking about.


texts: Loading GeoTiff with specified color band from ArcGIS to QGIS?

texts:
I have just started experimenting with QGIS. I have a GeoTIFF that I loaded in QGIS and used Singleband psuedocolor to create nice color codings for the map. However, my employer has specific colorband stored in his computer that he wants me to use. I don't know how can I export his specific colorband from ArcGIS and use it in my QGIS.

So far, I tried using the .clr file, but it didn't help me out.


texts: Problem with running and installation of QGIS

texts: Yesterday I was trying to upload a map file .tar in order to learn how to perform land cover classification with the semi-automatic plugin created by Luca Congedo. As I Woke up I saw that the QGIS application tab wasn't opening so I closed it and try to launch it again. The application launched but I couldn't open the tab so I could work on it. After restarting my PC I decided to uninstall and install it again. Recently I'm trying to download but first I get error messages so I canceled the installation and try to start all over again, it happened again so now I paused it and started again it did work but now I´m stuck at 0B/s speed.


texts: Calculating distance to the closest feature from different layer using Fiel Calculator

texts:
I am looking for the distance of the closest point (Layer 1) to my atlas object (Layer 2). I want to use the expression builder.
Example: I have a House, which is also the atlas_feature, on 'Layer A'.
I have bus stations which are on 'Layer B'. I want to get the distances to the closest Bus Station on 'Layer B'.
For getting the Name of the closest Busstation on 'Layer B' I am using reffunctions geomnearest. However, I also need the distance to the closest one.


texts: Can we insert arbitrary HTML layers the ESRI ArcGIS JavaScript API?

texts: We would like to implement a completely custom visualization to be included in a map, whose other layers we'd like the JavaScript API to control. This custom visualization would be updated by hand in response to map events, and will be constructed by hand from HTML elements.  That is, it is not a layer with simple feature or image behavior as could be created by deriving a custom image, feature, or tile layer. Other than reverse-engineering the DOM constructed by the JavaScript API, is there any support for this cooperation in the API?


texts: Is GeoNetwork an appropriate choice when I need to automatically create large amounts of metadata?

texts:
I have inherited a GeoServer installation publishes ~13,000 - 18,000 WMS layers in a dozen or so datastores on a 12 hour schedule.  Other than what's in the GetCapabilities document, there is currently no metadata.

There  is plenty of data scattered about, most of it is already in a format that is pretty machine readable, so the thought was to automatically create matching metadata for each layer published.

The question is: is GeoNetwork an effective platform for automating the storage & update of all this metadata, or should I:

write fgdc / iso xml files to the file system
do something else


texts: ArcGIS Online secure service access- user roles

texts:
I've added a secure map service to my organization's AGOL account. I was prompted to enter a username and password while adding the service, and did so, choosing to not require a log in each time. From this service, I've created a WebApp. I've added the App to a group.

As the administrator, I can successfully view this content on the WebApp. However, users in the group that do not have their roles set as "Administrator" cannot view this content. The basemap displays, but the content is nonexistent. If I, as the administrator, change a user's role to "Administrator", then the content is available to them. Any other role, however, does not display the content.

For good reason, I do not want to give a user Admin privileges. What role other than Admin will allow a user to view secure content? Or, is there a custom role I can create where I can allow access to secure services without granting Admin privileges?


texts:
I am new to the GIS-World and just learned a couple of tools. Now I want to start with my actual project, but since i don't have a overall understanding I would have to hear your opinions about the feasability of my ideas.

I have a lot of data points which refer to car drives getting tracked by GPS. I want to group points which are really close to each other to be able to compare their sensor values.

First I thought of rounding their values to something around 5 decimal places to group them, but then the points would not lie on the street anymore.

My actual plan is to load german street data from OpenStreetMap-Data using osm2pgsql into PostGIS. Then I want to calculate the distance of my data points to the closest street using somehow (maybe ST_Distance). Then I would like to find some algorithm to group points on the street by some flexibel interval (e.g. 5 meters).

Do you think this might work? Any other tools or data I should have a look for? Just trying to minimize programming time invest into some system which might not be able to work.


texts: How to read the Microstation element attribute values from .Net C#?

texts:
I am using following code for read the element attribute from index.But it shown the error like Unable to cast object of type 'System.Object[]' to type 'MicroStationDGN.DataBlock[]'.Is there any problem in my code?

 private string readAttribute(Element ele,int attribIndex)
    {
        String getStr = "";
        DataBlock[] dtBck = (DataBlock[])ele.GetUserAttributeData(attribIndex);
        if (dtBck.Length != -1)
        {
            dtBck[0].CopyString(getStr, false);
        }
        return getStr;
    }


texts: QGIS georeferencer - GCP projection problem

texts:
In short, the problem is that when Geoereferencer is reopened with a previously edited raster, all the previously saved control points have moved off the image to a location thousands of kilometres away. (QGIS 2.6, Windows 64 bit).

In Georeferencer, I create a number of GCPs (i.e. control points) by referencing to objects in the map canvas. After adding points, I start the georeferencing calculation. Operation is successful, the raster appears on the map canvas, along with the control points. The control points are also visible in the Georeferencer window, still overlying the points chosen for georeferencing. I save the control points and close the Georeferencer. All good so far.

Now for the problem. Having closed the Georeferencer, I reopen it, and open the new, georeferenced, raster image. But now, the control points in the Georeferencer window no longer overlay the raster in the Georeferencer window. However, they do appear correctly on the map canvas (the main QGIS screen).  In the Georeferencer window, they have in fact displayed, but half a world away to the south. That is, if the raster covers a rectangle with coordinates 300,000, 5834300 ; 362100, 5782700, the control points occupy an area bounded by coordinates 0,0 ; 6800, -5900. These small coordinates seem to correspond to the number of pixels in the raster image. So, it is no longer practically possible to edit the GCPs.

How do I solve the problem?

See screenshots below.

Image 1: Georeferencer window, and map canvas. After reopening Georeferencer, GCPs do not show in window, but show on canvas (the small red dots).

Image 2: Zooming out in Georeferencer window shows the GCPs are there but far away.


texts: Spatial Join also Based on Attributes using QGIS

texts:
I have two data sets:

points with a timestamp
polygons representing neighbourhoods

I need to compute how many points are in each neighbourhood polygon, based on a time range. So in other words, how many points are in zone A between 9 am and 10 am and then how many are there between 10 am and 11 am.
Currently, I am separating the points into separate layers - one layer for 9 am to 10 am then I do a spatial join. Then another layer for 10 am - 11 am and another join. This is exceptionally tedious and insufficient to do this same analysis across a week for example.
Any suggestions?


texts: Difference between types of Vegetation indices

texts:
I am currently comparing different types of Vegetation Indices, in this example I am using a NDVI formula and a TCARI formula. The NDVI result came out as expected with values ranging from +1 to -1. However the TCARI result came out with a different range. Kindly refer to the image below, illustrating the results of the TCARI Index before and after a bergwind hit a citrus farm. These images were obtained through use of an UAV and Micasense Rededge multispectral camera.

My question is, how do you interpret the values presented by TCARI? and how do you compare those values to that of NDVI?

I know the common range for green vegetation of NDVI is between 0.2 and 0.8, higher being more favourable. But how do I interpret TCARI without knowing the common range of values indicative of healthy vegetation?

.


texts: How to calculate final image size from wms?

texts:
I would like to retrieve an image from a WMS server (see the link below) with python OWSLib (or any other methods).

I know that resolution is 1m per pixel.

For example, I want to get bbox=(54.2652, 18.5889, 54.4189, 18.7866) (epsg 4326).

How can I calculate image size in pixel?

http://mapy.geoportal.gov.pl:80/wss/service/img/guest/ORTO/MapServer/WMSServer


texts: gdal CoordinateTransformation.TransformPoint result is not same gdal.Warp

texts:
I produce NDVI data from Landsat8 (UTM). End of the calculation, I want to export result EPSG:4326 so I used code that is below. I realized my new coordinates and qgis gdal warp result are different. lrx and lry are coordinates calculated wrong. Any idea?

g = gdal.Open(input_path)
red = g.ReadAsArray()
s_srs = g.GetProjectionRef()
osng = osr.SpatialReference ()
osng.SetFromUserInput ( s_srs )
wgs84 = osr.SpatialReference ()
wgs84.ImportFromEPSG ( 4326 )
tx = osr.CoordinateTransformation(osng,wgs84)
geo_t = g.GetGeoTransform ()
x_size = g.RasterXSize # Raster xsize
y_size = g.RasterYSize
(ulx, uly, ulz ) = tx.TransformPoint( geo_t[0], geo_t[3])
x=geo_t[0] + geo_t[1]*x_size
y=geo_t[3] + geo_t[5]*y_size
(lrx, lry, lrz ) = tx.TransformPoint(x,y)

#gdal TransformPoint result
38.839209332307 38.528922768342554 0.0
41.406489023769865 36.39800679336268 0.0
#qgis gdalwarp result
#Extent 38.8392093323070071,36.3980038478590018 :
41.4749801893845103,38.5290331355726110


texts: How to use GRASS Generalize to reduce GADM data size to 1/5 preserving topology?

texts:
The GADM is a great resource from global boundaries of countries, states and municipalities. The amount of detail is suitable for many analyses but when it comes to mapping or displaying it is often sufficient to display a simplified version.

I am struggling to create a simplified set of the GADM level 0 boundaries in QGIS using the grass generalize (v.generalize) function.

Objective:

Reduce total size of the file from ~500MB to <100MB
Create matching boundaries between countries. No weird empty space between polygons

so far I haven't been able to find the right settings of the function.

I am aware of the following resources:

tutorial
function page
Blog describing my probem


texts: ArcGIS Labeling Expression when truncation needed to fit in some polygons?

texts:
What would be the expression coding for if the address is 14802. But if the whole address cannot fit inside the building dimessions, then it should be truncated to the last two digits - 02. I am trying to do this ...

As seen in my map, I am trying to take down the congestion of the map with the numbers and trying to still be able to label each building for a firedepartment map.

Any ideas?


texts: What does enable sync do when uploading feature layer from ArcGIS Pro to ArcGIS Online?

texts:
Whenever I upload a feature layer from ArcGIS Pro to ArcGIS Online (AGOL), I see this option:

I want to know what this option means.

Does it mean there is a live connection between the feature layer in Pro and the uploaded layher in AGOL?


texts: Can I split my data based on the polygons in a different layer?

texts: Carto has a built-in dataset with polygon data for UK administrative regions. I have put another layer over this, containing several thousand data points (postcodes). Is there an easy way to find out how many post-code points fell into each administrative region, or do the layers not interact in this way?


texts: Iterating through a Map Series Layer using two map frames with ArcPy for ArcGIS Pro

texts:
I have one layout with 2 map frames on it. I want to create a map series where when I export to PDF, frame 1 iterates the odd numbered pages and frame 2 iterates the even numbered pages. The problem I have is at line 48 where I need to iterate on the second map frame. I can't figure out how to do it. Here is my code:
# Purpose: This script is designed to show 2 pages from a Map Series on
#          one Layout.  The first map frame runs through the odd numbered
#          pages in the series, while the second frame runs through the even
#          numbered pages.It then combines all pages into a single, multipage PDF.

import arcpy, os

#Filepath variables are machine specific. May need to change to execute sample
outputFolder = 'C:\\'   ###CHANGE if needed
pdfFileName = 'InsertNameHere.pdf'              ###CHANGE if needed

#Reference the current project, MUST be run in the ArcGIS Pro application
p = arcpy.mp.ArcGISProject('current')

#Reference the appropriate Map frames, Maps, Layers, and Layout
m = p.listMaps('MapNameHere')[0]                   ###CHANGE if needed
lyt = p.listLayouts("MapLayoutNameHere")[0]       ###CHANGE if needed
df1 = lyt.listElements("MAPFRAME_ELEMENT","MapFrameNameHere 01")[0]  ###CHANGE if needed
df2 = lyt.listElements("MAPFRAME_ELEMENT","MapFrameNameHere 02")[0]  ###CHANGE if needed
lyr = m.listLayers("IndexLayerNameHere")[0]   ###CHANGE if needed

####################

features = {row[0]:[row[1], row[2], row[3]] for row in arcpy.da.SearchCursor
       (lyr,['Id', 'SHAPE@','Scale','Name'])}

#####################

#References the active map series
ms = lyt.mapSeries

#Create a new PDF that pages will be added into
#First remove the file if it already exists
if os.path.exists(os.path.join(outputFolder, pdfFileName)):
  os.remove(os.path.join(outputFolder, pdfFileName))
pdf = arcpy.mp.PDFDocumentCreate(os.path.join(outputFolder, time.strftime("%Y%m%d") + pdfFileName))

#Iterate through the MapSeries
if ms.enabled:
    #This part itertates on the first map frame (Odd Numbers only)
    for pageNum in range(1,ms.pageCount,2):
        ms.currentPageNumber = pageNum

        #This part is supposed to iterate through the maps series on the second map frame (Even Numbers Only)
        arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",features[pageNum+1])
        df2.camera.setExtent(df2.getLayerExtent(lyr, False, True))

        #Export each page and append into the new PDF document
        print('Exporting and appending: ' + pageName)
        lyt.exportToPDF(os.path.join(outputFolder, pageName +  '.pdf'))
        pdf.appendPages(os.path.join(outputFolder, pageName +  '.pdf'))

        #Remove the single, exported page after appended into final PDF
        os.remove(os.path.join(outputFolder, pageName + '.pdf'))

#Complete the creation of the PDF and commit to file
pdf.saveAndClose()

Can someone help me figure out what I'm doing wrong?


texts: Copy columns from one layer to another QGIS 3.x

texts: I couldn't find any suitable answer for my struggle. I have two layers with different objects (buildings from different areas). I want them to have same columns, but creating one by one is very time consuming, so I was wondering if there's easier solution to copy empty columns from one layer to another?


texts: How to apply style in GeoServer REST API

texts:
My goal is to updade style through the REST API.

My question is, how to apply style through REST API.

When I am updating existed style via rest API, it only updates content.
It needs to be applied for loading style on layer.
My code below:

This command is just updating the content of style but does not apply to layer

 curl -u admin:geoserver -XPUT -H "Content-type: 'application/vnd.geoserver.geocss+css" -d @data_clone.css http://localhost:8080/geoserver/rest/workspaces/sde/styles/style_clone.css

It updates CSS fie but does not update SLD file


texts: Selection of points not showing on top (QGIS 2.2)

texts:
is there a way to bring the display of the selected objects on top? I have a large numer of points with a lot of them on top of eachother. When I make a selection through an expression and zoom to the selected points there is nothing to see. I have not found any option to bring the selected elements all the way to the top.

Setting the point symbolgy to "exploding" display does not help here either.


texts: Great circles in QGIS and export in 3857 web map

texts:
My plan is to visualize great circles in a Leaflet web map.
So coming from this great postgis article and R version  I thought it would work like this in QGIS:

create a line from two points,
store polyline in EPSG:4326
add many vertexes using the densify algorithm
store the densified line as 4326
put it into Leaflet and you should have nice "curves" on a map

what comes out is this:

And here a comparison of densified and normal lines in QGIS where you can see some small curvature of lines:

Is it just a problem with the EPSG:3857 or am I interpreting the post from @underdark wrong?


texts: OpenLayers 3, hover (select interaction with pointMove condition) on overlapping features

texts:
It is not recommended to have multiple select interactions in Openlayers. I can create a hover interaction using the select function with the pointMove condition to highlight features and display a popup if required. The interaction targets multiple layers and the styling is defined by the source of the feature. The problem which I am having is with overlapping features. In particular with lines and points from separate layers which intersect polygon features.

If for example I run with cursor along a line into a polygon the line stays highlighted as expected. The same is true if I highlight a polygon but then want to highlight a point inside the polygon or a line which crosses the polygon. This seems to be impossible with the current library. The pointer condition is still true for the polygon geometry but I would like to assign an index to the other overlays with point or line features to highlight these features if the pointerMove condition validates.

Here is jsfiddle to show the issue.


texts: Calculating angle of polyline at midpoint using Python/ArcPy?

texts:
Trying to calculate the angle of a line at its mid point using python; this will be used to display flow direction on a water network.

I've developed this script, which works in terms of generating the correct output, however, the angle returned is incorrect. I am currenlty getting some negative values, where I would like to have outputs of 0 to 360, in either an Arithmetic or Geographic format.

#OperationAlligator.py
import arcpy, math, datetime, numpy
from arcpy import env

print ("starting")
start = datetime.datetime.now() # for calculating time of process

#setting the test sources
flowArrow = r"J:\PYTHON\Flow_Direction.gdb\test" # pre-added fields "angle"
pipes = r"J:\PYTHON\Flow_Direction.gdb\Pipes"

#UpdateCursor
rows = arcpy.UpdateCursor(pipes)
updateRows_FlowArrow = arcpy.UpdateCursor(flowArrow)

#housekeeping
arcpy.DeleteFeatures_management(flowArrow)
arcpy.env.overwriteOutput = True

#Creating the physical xy locations in memory
shapeName = arcpy.Describe(pipes).shapeFieldName
for row in rows:
    feat = row.getValue(shapeName)
    x = feat.positionAlongLine(0.51,True).firstPoint.X - feat.positionAlongLine(0.49,True).firstPoint.X
    y = feat.positionAlongLine(0.51,True).firstPoint.Y - feat.positionAlongLine(0.49,True).firstPoint.Y
    radian = math.atan2(y,x)

    degrees = radian * 180 / math.pi
    print degrees
    row.angle = degrees
    rows.updateRow(row)

#Generating the mid point flow arrow
    out_cursor = arcpy.InsertCursor(flowArrow)
    midpoint = feat.positionAlongLine(0.50,True).firstPoint
    outRow = out_cursor.newRow()
    outRow.Shape = midpoint
    outRow.angle = row.angle
    out_cursor.insertRow(outRow)

del rows, row, updateRows_FlowArrow, outRow, out_cursor
print "Done in ",datetime.datetime.now() - start, " seconds"


texts: Time series export from ArcGIS Desktop to KML not working in Google Earth?

texts:
I am having trouble with the KML export. My time series will not work properly in google earth, it has time start in 0001 when it should start in the year 2001.

My incidents time is coded YYYYMMDDHHSS do you think if I cut it back to YYYYMMDD that it would work?

It is a geocoded series of terrorist attacks in Afghanistan.

Works fine in arcmap 10 and arcglobe 10.


texts: GRASS and SAGA function not working with QGIS 2.18.12

texts:
I was using QGIS 2.18.7 version where SAGA and GRASS was working fine. On updating to 2.18.10 and now 2.18.12, both SAGA and GRASS has stopped functioning.
QGIS is installed on win 7 64 bit system.

Does anyone have a solution to this?


texts: Getting a constant background color for transparency

texts:
I am creating a series of tif files at 1:10k scale over a vast area (760 tiles) using arcpy.mapping.ExportToTIFF and then stitch them together to create a single ecw. The output looks ok but when I try to turn on transparency I find that the background "white" is not pure white (255,255,255) but values that "look white" like "254,254,252" or "254, 251, 250)

Is there anyway of setting a tolerance in the export so that any value between say 250 and 255) get re-described as 255,255,255? or any no-data areas get 255,255,255?

Alternatively would putting a background color (like green screening) that's not on the mxd work better than keeping the background undefined?

I am posting this and also researching the above options and others as they come to me as this is urgent.

EDIT 1
When I try jpg2000 I get better results but some white's are not correct. I am running the following command to create the tif

arcpy.mapping.ExportToTIFF(mxd,r"tests\SF\100k\SF100k_" + str(ddp.pageRow.PageName) + ".tif",df,9306,9306*ar,1800,True)

The issue is with the generated tif's as the "white" changes between tiles so when I stitch them in Global Mapper the "white" is different on each file. Would a raster calculator operation be able to clean this up? or can I fix it in the generation of the tif (ideally)? I need to run a similar process on 16 other map series at a few different scales so we are looking at over 10,000 individual tif tiles in 16 groups so the latter is preferred.

EDIT 2:

PNG image

Same area TIF image

Any idea on why they look so different when it's exactly the same data? It's interesting that the symbology for the PNG is based on the data view while the TIF is based on the layout view. The PNG is much closer to the intended symbology.

Can Geoserver deal with Transparency in a different way to ArcGIS where the above issues then become irrelevant.

Thanks all,


texts: Creating buffers around locations when using a Geographic Coordinate System

texts:
Background

I am trying to make students familiarize with some some basic GIS procedures (using QGIS 2.18), like creating simple maps featuring points that represent the location of shops in a given area. I am making them use OSM as base map, which (if I am not mistaken) uses a Geographic Coordinate System (wgs1984).

Problem

I anticipate a problem when it will come to create buffers around those locations. The problem I envisage is related to the coordinate system used by the base map against which the point shapefile has been created.

Question

How to meaningfully create buffers when using a GCS? Is there a way to reproject the basemap to a projected coordinate system so that buffers will meaningfully employ metric units?


texts: arcmap - automatically change name of map

texts:
I would like to insert text place in arcmap with street name. For every I have to plot one map. I have hundreds of them in shp file and it would be great to avoid one by one manual typing street name.
Attribute table contain name for particular street i.e.:

I do not want to place one by one name of the street, then plot, move to next street area, change name, plot etc.
Is it possible to autoamatically change name of street in txt when I have different area in Layout Viev?


texts: How to create shape in Shapely from an SVG  element?

texts: AFAIK, you need to provide a list of coordinates to Shapely in order to create a shape but that's not how SVG <path> element works. Is there a way to create a Shapely shape using an SVG <path> element?


texts: Export raster over a loop in Google Earth Engine

texts:
I would like to export raster in a foreach loop with :

var landCover = ee.Image('MODIS/006/MCD12Q1/2016_01_01').select('LC_Prop1');
var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');

var Country_array = ['Kenya', 'Senegal',
                     'Mozambique', 'Comoros',
                     'France', 'Cambodia',
                     'Colombia','Gabon'];
Country_array.forEach(function(c) {
  print(c);
  var sel = countries.filter(ee.Filter.eq('country_na', c));
    var modis2017 = ee.Image(landCover)
                .clip(sel);

    Export.image.toDrive({
      image : modis2017,
      description: 'image'.concat('_', c),
      //maxPixels: 3784216672400,
      scale: 1000,
      crs : 'EPSG:4326',
      fileFormat: 'GeoTIFF'
      });

   // display
  Map.addLayer(sel, {}, 'shp'.concat('_', c));
  Map.addLayer(modis2017, {}, 'r'.concat('_', c));

});

but it produces empty (black) images. I don't understand why.


texts: ArcGIS Pro error 000732: "Input Layer: Dataset xyz does not exist or is not supported" when using arcpy.SaveToLayerFile_management

texts:
Below is the code I wrote in ArcPy in ArcGIS Pro. The layer from "arcpy.MakeFeatureLayer_management" is created (name changes in each iteration as "Layer_%s" %(FIPS)). But "arcpy.SaveToLayerFile_management" is not accepting that created layer and giving the error-000732. Ultimately I want to apply symbology from one layer to many layers generated from "SaveToLayerFile". Is there a work around it?

import arcpy

arcpy.env.addOutputsToMap = 0
arcpy.env.overwriteOutput = True

arcpy.env.workspace = "C:/Documents/Data/ArcGIS_WorkFiles/Metric_1/Monitors_Joined_With_PM_2000Size_AllCounties"
PM_Shapefiles = arcpy.ListFeatureClasses()

in_symbology_layer = "C:/Documents/Data/ArcGIS_WorkFiles/Metric_1/Monitors_Joined_With_PM_2000Size_AllCounties/Source_Symbology_Layer/Size2000_MaxTotalPM25_Layer.lyrx"

for shp in PM_Shapefiles:
    FIPS = shp.replace("PM25Monitors_JoinedInMAPPM_2000Size_", "")
    FIPS = FIPS.replace(".shp","")
    FIPS = int(FIPS)
    print ("FIPS=%s" %(FIPS))

    arcpy.Delete_management("Layer_%s" %(FIPS))
    where_clause = '"TotalPM25">0.0001'
    arcpy.MakeFeatureLayer_management("C:/Users/mthind/Documents/EMRE_Project/Data/ArcGIS_WorkFiles/Metric_1/Monitors_Joined_With_PM_2000Size_AllCounties/%s" %(shp), "Layer_%s" %(FIPS), where_clause)
    arcpy.SaveToLayerFile_management("Layer_%s" %(FIPS), "Layer_Files_2000Size/Symbology_%s.lyrx" %(FIPS))

    symbologyFields = [["VALUE_FIELD", "PM_MTInMAP", "PM_MTInMAP"],
                   ["COLOR_EXPRESSION_FIELD", "TotalPM25", "TotalPM25"]]

    arcpy.ApplySymbologyFromLayer_management("Layer_Files_2000Size/Symbology_%s.lyrx" %(FIPS), in_symbology_layer, symbologyFields)

Edit code to response @ 8:31 ET

import arcpy
import os

arcpy.env.addOutputsToMap = 0
arcpy.env.overwriteOutput = True

arcpy.env.workspace = "C:/Users/mthind/Documents/EMRE_Project/Data/ArcGIS_WorkFiles/Metric_1/Monitors_Joined_With_PM_2000Size_AllCounties"
PM_Shapefiles = arcpy.ListFeatureClasses()

in_symbology_layer = "C:/Users/mthind/Documents/EMRE_Project/Data/ArcGIS_WorkFiles/Metric_1/Monitors_Joined_With_PM_2000Size_AllCounties/Source_Symbology_Layer/Size2000_MaxTotalPM25_Layer.lyrx"

#workspace = "C:/Users/mthind/Documents/EMRE Project/Data/ArcGIS_WorkFiles/Metric_1/Monitors_Joined_With_PM_2000Size_AllCounties/OK"

for shp in PM_Shapefiles:
    FIPS = shp.replace("PM25Monitors_JoinedInMAPPM_2000Size_", "")
    FIPS = FIPS.replace(".shp","")
    FIPS = int(FIPS)
    print ("FIPS=%s" %(FIPS))

    arcpy.Delete_management("Layer_%s" %(FIPS))
    where_clause = '"TotalPM25">0.0001'
    arcpy.MakeFeatureLayer_management("C:/Users/mthind/Documents/EMRE_Project/Data/ArcGIS_WorkFiles/Metric_1/Monitors_Joined_With_PM_2000Size_AllCounties/%s" %(shp), "Layer_%s" %(FIPS), where_clause)

    if arcpy.Exists("Layer_%s" %(FIPS)):
        print("layer exists")
    else:
        print("Does not exist")

    if os.path.exists("Layer_Files_2000Size/Symbology_%s.lyrx" %(FIPS)):
        arcpy.Delete_management("Layer_Files_2000Size/Symbology_%s.lyrx" %(FIPS))

    arcpy.SaveToLayerFile_management("Layer_%s" %(FIPS), "Layer_Files_2000Size/Symbology_%s.lyrx" %(FIPS))

    symbologyFields = [["VALUE_FIELD", "PM_MTInMAP", "PM_MTInMAP"],
                   ["COLOR_EXPRESSION_FIELD", "TotalPM25", "TotalPM25"]]

    arcpy.ApplySymbologyFromLayer_management("Layer_Files_2000Size/Symbology_%s.lyrx" %(FIPS), in_symbology_layer, symbologyFields)


texts: Export georeferenced image overlay from KML?

texts:
How might one export image(s) from image overlays in a KML file to a standard geo-referenced raster format like geotiff?

Use case: NASA's Earth Observatory has a time series of the changing Athabasca Oil Sands, (kml). It's hard to get a personally meaningful sense of scale from the line scale bar on those images, or from viewing the scale bar in Google Earth. So I'd like to superimpose things I'm familiar with, but to do that I need to have them in a more manipulable format, with scale.

(I know that for this particular use case I can get the source geo-referenced images from NASA and just sidestep the whole conversion issue. However it's the general question I'm interested in answers for.)


texts: Source for all census tables available at each census geography/summary level

texts:
I'm familiar with using factfinder2 and generally know my way around census products. But one thing that I've never found is a document that has all the tables from major census products like the dicennial census or ACS and the level of census geography at which it's available.

Anyone ever come across such a thing?


texts: Suburb shapefiles from the Australian Bureau of Statistics (ABS) hass missing suburbs, unsure how to fill the holes

texts:
I downloaded the latest suburb shapefiles from the latest census at

http://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/1270.0.55.003July%202011

from the link

"State Suburbs ASGS Non ABS Structures Ed 2011 Digital Boundaries in ESRI Shapefile Format".

It has some missing polygons as when I open the shapefile in ArcGis for example it has holes around the Canberra region.

Im unsure how to best fill these holes are if there are tools out there to help me populate the holes with the correct polygons.

Thanks


texts: Can't upload Geotagged Photo folder in QGIS via Photo2KMZ?

texts:
I can't upload geotagged photos folder in QGIS via photo2KMZ. It's showing
error as follows

An error has occurred while executing Python code:

KeyError: 'DateTime'
Traceback (most recent call last):
  File "/home/yash/.qgis2/python/plugins/photo_to_kmz/photo_to_kmz.py", line 173, in run
    dt1,dt2=a['DateTime'].split()
KeyError: 'DateTime'

Python version: 2.7.6 (default, Oct 26 2016, 20:33:43) [GCC 4.8.4]
QGIS version: 2.14.12-Essen Essen, exported

Python Path:
/usr/share/qgis/python/plugins/processing
/usr/share/qgis/python
/home/yash/.qgis2/python
/home/yash/.qgis2/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python2.7
/usr/lib/python2.7/plat-x86_64-linux-gnu
/usr/lib/python2.7/lib-tk
/usr/lib/python2.7/lib-old
/usr/lib/python2.7/lib-dynload
/usr/local/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages/PILcompat
/usr/lib/python2.7/dist-packages/gtk-2.0
/usr/lib/pymodules/python2.7
/usr/lib/python2.7/dist-packages/ubuntu-sso-client
/home/yash/.qgis2//python
/usr/share/qgis/python/plugins/fTools/tools


texts:
I have created two checkboxes. If I check one it will zoom to a point with marker. If I uncheck the checkbox the marker removed. Similarly second one is working. But I want to do if the two checkboxes are checked at a time it will zoom at a position where both are visible. Again when one of it checkbox is unchecked it will zoom to the points best fit position. I know that fitbounds can solve this but couldn't imliment it. Here is my code I used. Any suggestion will be a great help. Thanks

<input type="checkbox" name="One" id="one" >One</input>

<input type="checkbox" name="Two" id="two">Two</input>

 $('#one').click(function(){
     if($("#one").prop('checked') === true) {
         var marker = L.marker([32.304255,-64.752098]).addTo(map).bindPopup('One');
         map.setView(new L.LatLng(32.304255,-64.752098), 16);
         marker.on('dblclick', function(){
         marker.openPopup();
         });
     }
 });
 $('#two').click(function(){
     if (this.checked) {
         var marker= L.marker([32.319343,-64.728066]).addTo(map).bindPopup('Two');
         map.setView(new L.LatLng(32.319343,-64.728066), 16);
         marker.on('dblclick', function(){
         marker.openPopup();
         });
      }
 });

var xLayer =[];
function clearLayer(){
for(var i =0; i < xLayer.length; i++)
    map.removeLayer(xLayer[i]);
    xLayer = [];
};

$('#one').click(function(){
    if($("#one").prop('checked') === false) {
    clearLayer();
    }
});
$('#two').click(function(){
    if($("#two").prop('checked') === false) {
    clearLayer();
    }
});


texts: esri/dijit/editing/AttachmentEditor use proxy when viewing attachments

texts:
I'm trying to figure out why the ArcGIS JavaScript API AttachmentEditor does not use a proxy when viewing attachments.  I have my proxy set up correctly, and it uses it in every other situation but this.

I'm coming up empty on googling for the answer and have found these threads on Esri which don't give any resolution to it.

https://community.esri.com/thread/47357
https://community.esri.com/thread/176862-attachmenteditor-with-proxy

Has anyone had experience with getting the view attachments to go through a proxy?


texts: Calculating Proportion of Land Cover Classes with moving window around point in R?

texts:
I am currently using a raster that is a land cover class raster (with 50 different land cover classes at a resolution of 30m). I would like to obtain 50 different rasters (essentially predictors for a future regression analysis) that are data of a continuous nature.

Previously, I used the focal statistics tool in ArcMap but I ended up obtaining either a 0 or a 1 based on areas where a particular land cover type was present.

See image below (All areas of White are 1 and Areas in black are 0):

Is there a way to obtain multiple rasters from a single land cover class raster, with the proportion of each land cover class calculated around a window of 5*5 from a particular point (instead of a 0 or a 1) in R?

I noticed that there is a focal function in the raster package.


texts: Clipping a raster dataset with different projections

texts:
I have downloaded the Corine Land Cover map 2018 (raster 100m) from Copernicus Services, which covers the European Continent.
As I only need to work within a certain area, I need to clip the original data from another dataset that I have, to work just with my AOI.
The work that I will need to do is about spatial and statistical analysis and modelling.

My question is concerning projections when clipping the data from the original dataset, and the resulting output rasters.

The original dataset projection is ETRS 1989 LAEA (EPSG 3035).
My AOI projection is ETRS 1989 Transverse Mercator. If its correct to say, I may add that the latest is the one that should be appropriate because is local. The issue comes next:

a) If keeping the original projection while clipping, the output dataset will have the pixels in a diagonal orientation, which is not what I want, and the Sum of the cell count will be 8,925,616.

b) If transforming the original projection into the local projection that is used in the region, the output dataset will have the pixels in the vertical orientation, which is what I want, but the sum of the cell count will be higher: 8,926,281

What is the best approach for the following scenario and why?


texts:
I am building a statistical analysis project using ASP.NET, Geoserver, and OpenLayers.
I need to let the user draw a polygon on OpenLayers then find geometries(villages etc.) that intersects this polygon.  Then get their attributes like temperature, population etc., make some mathematical calculations and present this information to user.

My question is can I do the following operation:

(SELECT population,temperature,humidity,height FROM villages WHERE geometry INTERSECTS userDrawnPolygon,
SELECT data1,data2,data3 FROM tableX WHERE geometry INTERSECTS userDrawnPolygon .... goes on for multiple tables)
with ECQL or CQL or some filter on Geoserver or do I need to develop this algorithm in my ASP.NET code - connecting to MSSQL?

Thanks already.


texts: Saving 'OUTPUT': from processing toolbox to a file path with unique file name

texts:
I'm using the QGIS Toolbox and Python Console to write a script for parsing data.

I want to do something like this:

output_path = open('C:/Users/EH/Desktop/PYQGIS_Model_Output') # Set location to save to
#Algorithm parameters
extractANN = processing.run('qgis:extractbyattribute', {
'FIELD': 'FlagANNtraining',
'INPUT': ColourCode['OUTPUT'],
'OPERATOR': 9, # Meaning extract by "FlagANNtraining" 'is not null'
'VALUE': '',
# Attempt to save output to path, with new file name
'OUTPUT': output_path + ANN.csv })

Below is a snipped of the ColourCode['OUTPUT'] - so the column is full of 1 or NULLs, whether it has been flagged or not.

The extract by attribute works fine when I'm saving as 'TEMPORARY_OUTPUT', so that isn't what I'm struggling with. I need to create a the kind of function where I can simply change the output path and have a ANN.csv for example saved there, and for each other 4 extracts it will be b.csv, c.csv for example.

EDIT: Thanks BERA for the snippet below.
Question #2

Code snippet that now works:

    infile = ColourCode['OUTPUT']
    output_path = ('C:/Users/EH/Desktop/PYQGIS_Model_Output')
    values = ['1']
    outprefix = 'flag_ANN_training_'

    for count, val in enumerate(values):
        extractANN = processing.run('qgis:extractbyattribute', {
        'FIELD': 'FlagANNtraining',
        'INPUT': infile,
        #'OPERATOR': 9,
        'VALUE': val,
        'OUTPUT': os.path.join(output_path, "{1}.csv".format(outprefix, count)) })

Why is the .csv that's created called '0', thought the outprefix would determine that?


texts: Exporting many small and spread out regions of interest of an image

texts:
I have a feature collection with more than 4500 polygons. Their summed area is ~100 km2, but they are spread out all over the USA. When I try to export a single image, I received an error message saying that the output exceeds the memory. I tried to export a smaller region that wouldn't exceed the limit and my output contains all the empty space between my polygons in black. I wonder if that could be the problem. I wouldn't mind having a separate file for each polygon but I don't know how to do that. Would really appreciate some advice on how to download this.

// Import Polygons
var polygons = ee.FeatureCollection("users/royjeanromain/walmart");

// Simplify Polygons
polygons = polygons.map(function(feature) {
  return feature.simplify(10);
});

// Get Data
var dataset = ee.ImageCollection("COPERNICUS/S2_SR")
  .filterDate('2019-10-01', '2019-12-31')               // date filter
  .filterBounds(polygons)
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10)); // no clouds

// Get the number of images.
var count = dataset.size();
print('Count: ', count);

// Reduce the collection with a median reducer.
var image = dataset.reduce(ee.Reducer.median());
print(image)

// Clip to polygon
var clipped_image = image.clipToCollection(polygons)
print(clipped_image);

// Only grab cetain bands
var reduced_median = clipped_image.select(['B4_median', 'B3_median','B2_median']);
print(reduced_median);

// Display a composite with a band combination chosen from:
// https://landsat.usgs.gov/how-do-landsat-8-band-combinations-differ-landsat-7-or-landsat-5-satellite-data
Map.addLayer(reduced_median, {bands: ['B4_median','B3_median','B2_median'], max: 10000},'also median');

//Export the image, specifying scale and region.
Export.image.toDrive({
  image: reduced_median,
  description: 'walmart',
  scale: 10,
  region: polygons,
});


texts: Update status of Bing Maps or Here Maps imagery for US area?

texts:
Does anyone know if any updates to Bing or Here (not sure what it is now) imagery is planned for the Southeastern Washington state area?

We are not highly populated but we have a 1/4-million urban area and we are growing quickly enough that some of our new developments still look like farm fields. I think current vintage shown in our area is 2011. I saw a good answer on the state of the data in Bing/Here Maps at BING Maps vs. HERE Maps, but not about the imagery. I ask because we use Mango Maps which defaults to Bing/Here as a basemap, and I am concerned about our web map presence looking "old" or 'dated'.


texts: Checking for raster coverage within polygon feature class using ArcPy?

texts:
I'm writing a Python toolbox and would like to build in a check for raster data within the user-supplied area of interest. This piece of code works most of the time:

if (str(arcpy.GetRasterProperties_management(chsaoi, "ALLNODATA")) == "1"):
    arcpy.AddWarning("No raster coverage in input OID {}.".format(row[1]))

But when the polygon is very large, this check seems to fail. Even when the entire input feature class has raster coverage, the "ALLNODATA" check returns '0'.

Can someone help me come up with a better way to do this?


texts: Datetime format in QGIS Points to Paths plugin

texts:
At 2.18.6, I am trying to use the PointstoPaths plugin.

My input layer has a datetime string field with the following format:

2015-01-25 13:45:50  [Please note the blank between date and time.]

However, the PointstoPaths window displays a date format example of:

%Y%m%dT%H:%M:%S    [Here, please note the 'T' between date and time]

Here's a screenshot:

If I use PointstoPaths with the following datetime format (which matches my data):

%Y%m%d %H:%M:%S  [with a blank!]

This error is returned:

If I subsequently edit my datetime string field to change the blank to a 'T', the tool executes perfectly.  Thus I am led to believe that the 'T' is required.  But I don't know what the 'T' represents.  A review of the Python datetime web page turned up nothing.

Sorry if this is long winded, but I don't understand the 'T'!


texts:
What is/are the best way(s) to host mbtiles?

I have tried:

hosting on my own server with tileserver.php: this works and is very straighforward and simple to set up, but doesn't seem to be super fast (maybe a problem with my server) and I fear that when the site opens to the world it will simply break down; (also, tileserver.php sends HTTP error 500 when a tile isn't found, which I find misleading?)
hosting on S3 (1 tile per S3 object): serving tiles is very efficient and fast, but uploading to S3 using tilecloud always raises an S3 error "connection reset by peer" in the middle of the process, so uploading never ends properly and more than half the tiles are not sent (and it seems that trying again just restarts everything from the beginning)

MapBox would be the obvious choice but their prices seem very high?

The mbtiles file is under 500 Mo for zoom levels [9,15], gets to 1 Go for levels [9,16] and 2.7 Go for levels [9,17]; I can live with [9,15].

S3 would be the simplest for me, if I can find a way to upload all tiles without error (or, with an error that can be debugged: S3 logging doesn't show anything for the error I'm experiencing).

Any other idea?


texts: Importing 3D models/wireframes into QGIS

texts:
I was wondering if Qgis2threejs is able to import 3D models (wireframes in 3D dxf) from GoCad, Datamine or any other 3D modelling software?

If not, is there any other plugins for QGIS that would be able to do this?


texts: Extend line (utility-water) to polygon (building)

texts: I have several thousand water lines that need to simply be extended to the edge of the nearest building. Some of my lines do touch or intersect the buildings and some don't. They are ranging from a few inches up to 50 feet from the edge. Is there a way to select the lines that aren't touching or intersecting a building and extend them so that they do touch or intersect a building?