texts: How to shift polygons by updating their details in R

texts:
I have two polygon shapefiles and I want to see to what extent the observed amount of overlap is due to chance. I'd like to randomly move the polygons in one shapefile to new locations, measure the total area of intersection, repeat this a thousand times or so to get a distribution of areas of intersection, and then measure how likely the observed area is compared to the distribution of random areas.

I'm stuck on moving the polygons: how can I shift polygons to a random new location in a bounding box? I'm using R, here's an attempt at a reproducible example:

require(sp)
require(splancs)
require(rgeos)
require(maptools)

# simple example polygons
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)

Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
Srs3 = Polygons(list(Sr3, Sr4), "s3/4")
SpP = SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3)
plot(SpP, col = 1:3, pbg="white")

# make bounding box to contain random relocations of polygons
bbx <- bboxx(bbox(SpP))
bbx <- rbind(bbx, bbx[1,])
bbx <- SpatialPolygons(list(Polygons(list(Polygon(bbx)), "bbx")))

# make random points in bounding box to use as centroids to sample from
rnd <- spsample(bbx, 100, type = "random")

# get centroids of polygons
cnt <- gCentroid(SpP, byid = TRUE)

Here's where I'm stuck...

# replace centroids of polygons with random points in bounding box
# in the hope that this will shift the polygons
out <- SpP # to store output
for(i in 1:length(SpP)) {
    # create a random offset to shift the polygon by
    offset <- unname(cnt[i,]@coords) - unname(rnd@coords[sample(1:length(rnd),1),])
    # shift the polygon by this random offset
    # this doesn't seem to be doing what I want...
    out[i,] <-  elide(SpP[i,], shift = offset)
}
plot(out) # not the desired result...

# check to see that something is happening...
sapply(1:length(SpP), function(k) SpP@polygons[[k]]@labpt)

The desired output is a shapefile with the polygons in new random locations in the bounding box. My general approach is based on this. What am I doing wrong and how can I fix it?


texts:
I have two raster with same resolution and dimension: water-occurrence index and rainfall (mm). I want to create a new raster of SLOPE (a) and INTERCEPT (b) by performing linear regression between every 3×3 pixels of the two rasters (rainfall and water-occurrence), such that each pixel of the SLOPE and INTERCEPT will hold the regression slope and intercept value obtained from linear regression of the corresponding 3×3 pixels in water-occurrence and rainfall that surround that pixel.
I adapt the code provide by @jeffrey-evans in this thread https://gis.stackexchange.com/a/278985
# Linear regression between every 3×3 pixels of two rasters
# Based on thread from SE: https://gis.stackexchange.com/questions/278979/linear-regression-between-every-3×3-pixels-between-two-rasters-using-r

library(raster)
library(gstat)
library(sp)

# Set working directory
setwd("/Users/mac/data/")

# Raster import and stacking
# Change the path variable
fs <- list.files(path="/Users/mac/data", pattern = ".tif$", full.names = TRUE)
r <- raster::stack(fs)

# Focal linear regression
w <- c(3,3) # 3x3 neighbour/group of cells

# Slope
lapse_1 <- r[[1]]
lapse_1[] <- NA
for( rl in 1:nrow(r) ) {
  v <- getValuesFocal(r[[1:2]], row=rl, nrows=1, ngb = w, array = FALSE)
  fit <- rep(NA,nrow(v[[1]]))
  for(i in 1:nrow(v[[1]]) ) {
    xy <- na.omit( data.frame(x=v[[1]][i,],y=v[[2]][i,]) )
    if( nrow(xy) > 4 ) {
      fit[i] <- coefficients(lm(as.numeric(xy$y) ~ as.numeric(xy$x)))[2]
      if(is.null(fit)) fit = 1
    } else {
      fit[i] <- NA
    }
  }
  lapse_1[rl,] <- fit
}

lapse_1
plot(lapse_1)

# Saving SLOPE into geotiff
writeRaster(lapse_1, filename=file.path("/Users/mac/output/Day1_01_Jan_Slope_3x3.tif"), format="GTiff", overwrite=TRUE)

# Intercept
lapse_2 <- r[[1]]
lapse_2[] <- NA
for( rl in 1:nrow(r) ) {
  v <- getValuesFocal(r[[1:2]], row=rl, nrows=1, ngb = w, array = FALSE)
  fit <- rep(NA,nrow(v[[1]]))
  for(i in 1:nrow(v[[1]]) ) {
    xy <- na.omit( data.frame(x=v[[1]][i,],y=v[[2]][i,]) )
    if( nrow(xy) > 4 ) {
      fit[i] <- coefficients(lm(as.numeric(xy$y) ~ as.numeric(xy$x)))[1]
      if(is.null(fit)) fit = 1
    } else {
      fit[i] <- NA
    }
  }
  lapse_2[rl,] <- fit
}
lapse_2
plot(lapse_2)

# Saving INTERCEPT into geotiff
writeRaster(lapse_2, filename=file.path("/Users/mac/output/Day1_01_Jan_Intercept_3x3.tif"), format="GTiff", overwrite=TRUE)

I am not familiar with R and programming, I am wondering whether above script is correct or not for calculating slope and intercept?


texts: Number of rows in free acount of cartoDB

texts:
How many rows in a table can have a free account?
and in the case of Magellan plan?


texts: Is it possible to save expression, so one doesn't have to retype?


texts: Edit xy domain of feature classes during import from boilerplate geodatabase schema (ArcGIS)

texts:
Situation:

I work for a company that creates geodatabases for a client that are restricted to a specific boilerplate schema (eg domains, subtypes, etc.).  During a recent draft review, the client complained that when zooming to full extent during an ArcMap session, the full extent was too far "out."

Desired solution:

Alter xy domain of empty feature classes (contained in feature datasets) to extent of interest area as they are imported into production geodatabase (I know that a feature dataset destination will ignore the xy domain environment setting).

Problem:

No geoprocessing tool (eg feature class-to feature class, import to geodatabase, project, define projection, etc.) will accomplish the task (inside or outside feature datasets).  The domain environment setting (advanced geodatabase) is ignored for all actions I've tried.


texts: Problem adding fields with Model Builder iteration

texts:
I'm trying to add 5 different fields to a bunch of shapefiles. I've made a model to automate it, which I'll post a picture of. It says it succeeded when it's finished running, but none of the shapefiles have the added fields in their dbf file. Any idea what I'm doing wrong?

Edit: The dbf files are in a subfolder of the folder I'm using for the iteration. I've made sure the recursive radio box is checked.


texts:
I'm new to Qgis. I have a script that runs successfully inside the Qgis editor, but it fails outside. It returns the alghelp, but fails on running the actual tool so it seems I have the external IDE communicating with modules properly. Ideas?

import sys
import os

from qgis.core import *

QgsApplication.setPrefixPath(r'C:\OSGeo4W64\bin', True)
app = QgsApplication([], True)

import processing
from processing.core.Processing import Processing
Processing.initialize()

alg = 'saga:fillsinksxxlwangliu'
processing.alghelp(alg)
srcDir = r'G:\stest'

for e in os.listdir(srcDir):
    if e.endswith('.tif'):
        print 'Processing ' + e
        processing.runalg(alg, srcDir + os.sep + e, 0.025, srcDir + os.sep + e[:-4] + 'filled.tif')

app.exit()
QgsApplication.exitQgis()

the traceback is

<type 'exceptions.AttributeError'>: 'NoneType' object has no attribute 'crs'


texts: ArcGIS API 3 Edit access ghost lines geometry

texts:
So my question is how does one get ghost line geometry while editing vertices using esri/toolbars/edit?
I would like to get the geometry of the ghost vertexes that are being drawn while moving the mouse in edit mode. However, it seems impossible for me to get it while in the mode of editing. Events on the map do not help (mouse-move or drag not recognized while editing) and the event on vertex-move gives me only information about the vertex geometry at its starting point the whole time while moving ghost lines(??).
Does anyone have any idea how could I get coordinates or graphic of ghost lines while editing?


texts: Does Geoserver support WMTS protocol? I find only WMS WFS and WFS-T support.


texts: Implementing INSPIRE Metadata with Geodatabase

texts:
I would like to begin by saying that I have never implemented INSPIRE metadata with my spatial datasets before.
I am undertaking a project which requires my geodatabase metadata to be in the standards format specified in Regulation 1205/2008/EC.
I believe I have most of the information required, for example 'Resource Title', 'Resource Abstract', 'Topic Category', 'Keyword', 'Geographic Bounding Box', 'Lineage', 'Access and use', Responsible organisation', 'Metadata point of contact', 'Metadata date', 'Metadata language', 'Coordinate reference system'.

My dataset are just vector points.
I am using ArcCatalog and ISO 19139 Metadata Implementation Specification.
I was wondering where I might put details of 'Conformity' and how should I go about recording 'Encoding' in the metadata.
Is ESRI software enough to create metadata which meets INSPIRE standards?


texts: Making geotiff and reprojecting to EPSG:4326 using GDAL?

texts:
My goal is to download a gif file from this webpage:

http://aviationweather.gov/data/iffdp/2104.gif

Then crop it to the known latlon bounds of 75N,45S,35W,70E, convert to geotif, and then project to EPSG:4326. The map says "Mer", but how does one figure out what EPSG to use in gdal when converting to geotiff? I've tried many: 3395,3857,3785, etc, but none seem to work properly. Here are the commands I've used:

gdal_translate -of GTiff -a_srs EPSG:3857 -a_ullr -35 75 70 -45 2104_cropped.png 2104_3857.tiff

gdalwarp -s_srs EPSG:3857 -t_srs EPSG:4326 2104_3857.tiff 2104_4326.tiff

I get no warning with gdalwarp, but the image never warps, it stays the same. When i overlay them on the maps that expect EPSG:4326, the maps to don't line up.


texts: arcpy.ArcSDESQLExecute with UPDATE and OUTPUT

texts:
Is there a way to get OUTPUT through arcpy.ArcSDESQLExecute? Whatever I do I get only True. When running UPDATE through MS SQL Server Management Studio I get OBJECTID of all updated features (rows):

UPDATE MY_TABLE SET MY_MINIMUM = 0 OUTPUT INSERTED.OBJECTID WHERE MY_MINIMUM < 5
DB result type: <type 'bool'>
DB result: True


texts: determine whether location of fish is wrong based on boat speed and time of last capture

texts:
I have a series of coordinates for fish caught from a boat at different datetimes and different trips. How do I determine whether the coordinates of a fish are likely to be incorrect (e.g. due to transcription error) based on time since last fish caught within that same trip and an assumed boat speed (say 10km/hour).

Here is a simple example dataset with 2 trips and two fish per trip.

library(sf)
library(ggplot2)

datetime <- ymd_hms('2017-05-13 14:00:00', tz = "Etc/GMT+8")
df <- data_frame(DateTimeCapture = c(datetime, datetime + minutes(35), datetime + days(2),
                                     datetime + days(2) + minutes(20)),
                 Trip = c('1', '1', '2', '2'),
                 Order = c(1, 2, 1, 2),
                 X = c(648635, 648700, 647778, 658889),
                 Y = c(5853151, 5853200, 5854292, 5870000))

# if you prefer to work in sf
df_sf <-  st_as_sf(df, coords = c('X', 'Y'), crs = 32610)

# quick plot
ggplot() +
  geom_point(data = df, aes(x = X, y = Y, color = Trip))

The distance between the two fish in the second trip is 19km:

st_distance(df_sf[3:4, ])
Units: m
         [,1]     [,2]
[1,]     0.00 19240.47
[2,] 19240.47     0.00

It is unlikely that a boat could travel 19km in 20 minutes. Thus this should be flagged as a possible error.

My preference is for solutions using sf, but may also accept solutions using sp. It has to be r-based solution.


texts: Adding button to tool dialog of Python Script tool?

texts:
I am getting stuck when I want to create a new python script in toolbox. The function I want to develop is to let user to input the user ID, then retrieve from the database the other values to allow user to do the updating.

In this case, I want add a "Retrieve" button in the script to allow user to get their details.

Is it possible?


texts: Efficient way to find geometries that intersect a given rectangle

texts:
It feels I'm missing something absolutely basic here OR there is something weird in the internals of Oracle's spatial indexing (I sincerely hope it is just me missing something basic or I am about to lose some faith in Oracle).

I'm trying to find which geometries intersect a given rectangle:

I have a rectangle, let's call it VIEWPORT with Latitude/Longitude 41.500, 12.500 and 42.000, 13.000

on-the-fly that is
SDO_GEOMETRY(2003,4326,null,SDO_ELEM_INFO_ARRAY(1,1003,3),mdsys.SDO_ORDINATE_ARRAY(12.500, 41.500, 13.000, 42.000))

Also I have in a table AREAS a bunch of geometries (polygons and multipolygons). Spatial indexed.

Now, I want to see which records from AREAS intersect that VIEWPORT. So, in my SELECT I should have a WHERE clause like this:

SDO_RELATE
(
    SDO_GEOMETRY(2003, 4326, null,SDO_ELEM_INFO_ARRAY(1,1003,3),mdsys.SDO_ORDINATE_ARRAY(12.500, 41.500, 13.000, 42.000)),
    "AREAS"."GEOMETRY",
    'mask=ANYINTERACT'
)='TRUE'

Well this can't execute because ORA-13226: interface not supported without a spatial index. That happens because the first parameter in SDO_RELATE is the on-the-fly SDO_GEOMETRY, which has no spatial index.

So let's flip the order of the arguments (which is not correct approach but anyway) just to mesaure the execution time. So the where clause goes like this:

SDO_RELATE
(
    "AREAS"."GEOMETRY",
    SDO_GEOMETRY(2003, 4326, null,SDO_ELEM_INFO_ARRAY(1,1003,3),mdsys.SDO_ORDINATE_ARRAY(12.500, 41.500, 13.000, 42.000)),
    'mask=ANYINTERACT'
)='TRUE'

"AREAS"."GEOMETRY" is spatial indexed, so no problem for SDO_RELATE to execute. I'm getting the correct results in 8.8 seconds. But, that's not efficient at all.

Now let's do something tricky. Instead of providing that VIEWPORT rectangle as a on-the-fly constructed SDO_GEOMETRY let's put it in a table, with spatial index. This way we can use it as the first argument on SDO_RELATE. So eventually we have:

WITH viewport AS
(
    SELECT * FROM "VIEWPORT_TEMP_TABLE" WHERE "ID" = 1
)
.
.
.
SDO_RELATE
(
    viewport."GEOMETRY",
    "AREAS"."GEOMETRY",
    'mask=ANYINTERACT'
)='TRUE'

This executes just in 3.8 sec, the difference is unexpectedly great!

So Oracle has the potential to perform this in 3.8 sec but only with this hacky way.

I think in this point, either I am ignoring something very basic and mis-using or mis-judging Oracle, or something is very "weird" in how Oracle performs intersections with a given SDO_GEOMETRY.

My question is: What is the correct/efficient way to find the geometries that intersect a given rectangle in Oracle?

As I demonstrated the fastest way is to create a temporary table, insert that given rectangle and rebuild its spatial index and then use it as first argument in the SDO_RELATE. This is way too hacky to be the correct way to perform something so simple like this.

So, what is the correct way to find geometries intersecting a given rectangle?

(In contrast, postgres has no problem with on-the-fly given rectangle in the very same scenario and that made me very sceptical about Oracle)

edit: running on

Oracle Database 11g Enterprise Edition Release 11.2.0.1.0 - 64bit Production
PL/SQL Release 11.2.0.1.0 - Production
"CORE   11.2.0.1.0  Production"
TNS for 64-bit Windows: Version 11.2.0.1.0 - Production
NLSRTL Version 11.2.0.1.0 - Production


texts: Converting raster pixels to vector cells?

texts:
Is there any way to convert raster pixels to vector cells?
I mean each pixel to a polygon.

I tried the SAGA tool 'Raster values to points' by choosing cells option instead of nodes option.

As a result, I got polygons but each one is not exactly  on the pixel. When I tried with nodes option, I got points but they are not in the center of the pixel.


texts: 'Spline' tool in ArcGIS Pro keeps dropping collinear points

texts:
I am trying to convert the WorldClim dataset from their original resolution of 30 arc-seconds (WGS1984) to 25m cells (OSGB British National Grid) for use in Maxent modelling.

I clipped the individual raster files (bioclim1, bioclim2 etc.) to my study area in south Wales and converted to OSGB using the 'Project raster' tool. I then converted the resulting files to point using 'Raster to point' and ran the 'Spline' tool with 25m as my new cell size.

The tool succeeded but included multiple warnings saying that 'collinear points have been dropped'. What does this mean? Is the tool designed to drop collinear points (which doesn't make much sense to me intellectually) or have I done something wrong?

For those with experience in Maxent modelling: Do I even need to go through all this or can I just use the original data in their original projection and resolution?


texts: Proj4 String for NAD83(2011) / Louisiana South (ftUS)

texts:
I am trying to convert NAD83(2011) / Louisiana South (ftUS) coordinates on survey maps to lat/lon using Proj4 in Python.

Could someone point me in the right direction?

I have verified the lat/lon on earthpoint.us site.  Their conversion seems to work.  This is my python code which contains the Proj4 string.

===============

from pyproj import Proj

#EPSG:6479 from epsg-registry.org
southla = '+proj=lcc +lat_1=30.7 +lat_2=29.3 +lat_0=28.5 +lon_0=-91.33333333333333 +x_0=3280833.3333 +y_0=0  +ellps=GRS80 +datum=NAD83 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs'

p=Proj(southla)

lonlat = p(3037416.890, 709211.549,inverse=True)

print "%s, %s" % (lonlat[1],lonlat[0])


texts: spatialite geometry vs PostGIS geography

texts:
When I worked with PostGIS I used geography type for geo objects with gps coordinates like latitude and longitude on the WGS-84 ellipsoid.
Now I have to use SQLite and libspatialite correspondingly.
I wonder, if I create column like this:

Select AddGeometryColumn ('my_table', 'Geometry', 4326, 'POINT', 2)

will queries with ST_Within/ST_Contains(geo_rect, my_table.geometry) find suitable
points if geo_rect covers north pole, or 0,180 longitude?

here how it works in postgis:

    geo_demo=# select ST_DWithin(ST_GeographyFromText('POLYGON((-179 53, 179 53, 179 50, -179 50, -179 53))'),  ST_GeographyFromText('POINT(-179 61)'), 1.0, true);
 st_dwithin
------------
 f
(1 row)

geo_demo=# select ST_DWithin(ST_GeographyFromText('POLYGON((-179 53, 179 53, 179 50, -179 50, -179 53))'),  ST_GeographyFromText('POINT(-179 51)'), 1.0, true);
 st_dwithin
------------
 t
(1 row)

geo_demo=#

as you see point with latitude 51,-179.5 inside geo box:

53, -179 <-> 53,179
50, -179 <->  50,179

but

spatialite> select Within(MakePoint(-179.5, 51, 4326), PolyFromText('POLYGON((-179 53, 179 53, 179 50, -179 50, -179 53))', 4326));
0

gives 0, but really point inside geo box.


texts:
I am relatively new to coding, but I am using fiona within python and would like to print the coordinates that make up the points of the polygon to later add to a separate list. However I am having a hard time printing just the coordinates, I do not wish to print the id, type, properties etc... My current code prints the points as well as the type, I would like it to print without the type. Is there a simple way to do this?

output_polygons = []
with fiona.open(input_shp) as polygons:
    for polygon in polygons:
        polygon['geometry']['coordinates'][0]
        output_polygons.append(polygon)
        pprint.pprint(polygon['geometry'])


texts: Python Error: "Data driven pages is not enabled", when they are enabled

texts:
I am exporting Data Driven Pages using arcpy (ArcMap 10.5, Python 2.7.13). I have used this script many times on this exact same MXD to export DDP to TIFF:

import arcpy

mxd = arcpy.mapping.MapDocument(r"C:\Users\juanzepeda\Desktop\WildCAD GIS\basemap.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):

  mxd.dataDrivenPages.currentPageID = pageNum

  arcpy.mapping.ExportToTIFF(mxd, r"C:\Users\juanzepeda\Desktop\WildCAD
  GIS\TIF10\basemap10.mxd" + str(pageNum) + ".tif", df, df_export_width=3000,
  df_export_height=3000, geoTIFF_tags=True)

del mxd

Once in a while, I get the error:

Runtime error
Traceback (most recent call last):
File "<string>", line 7, in <module>
File "c:\program files
(x86)\arcgis\desktop10.5\arcpy\arcpy\arcobjects\mixins.py", line 655, in
dataDrivenPages
return
convertArcObjectToPythonObject(self._arc_object.pageLayout.dataDrivenPages)
AttributeError: Data Driven Pages is not enabled on this map document.

I created a grid index of my features to use as the data driven layer. Then I make a selection of this layer (to break up large exports into chunks). Sometimes it works, sometimes I get that error. I try restarts, I have reset my %APPDATA% ESRI "user profile", as well as my default template in the %APPDATA% ArcMap 10.5 templates folder.

I try disabling DDP, then re enabling. I have tried Selection>Create layer from selected features, then using that layer for DDP, as well as exporting selected features to new shapefile, and using that for DDP. I don't know why this error occurs and how to fix it.

It happens when DDP is definitely enabled, you can click through them in Data View and it functions as expected. If I click on the attributes of DDP the correct layer and it's attributes show up. But somehow, the system isn't recognizing that DDP are enabled when using the ExportToTIFF.

I found this post on a different DDP export Attribute Error,but it is the only thing I can find evenly vaguely related to my issue. I did try removing the
    + ".tif" from my export file name like that post suggested but it does not help.

If I remove the layer that DDP is using, I get the warning that I am deleting the index layer being used for DDP. If I don't remove the layer and run my script, I get the error that says DDP are not enabled.


texts: Publish postgis raster based on column value in geoserver

texts:
I have a postgis raster table with the following columns:
rid, rast, year, plant

I created a wms layer in geoserver pointing to the table and successfully serve out the raster when the table only has one raster in it (one plant). The url used to do this looks like so:

http://<myserver>/geoserver/sde/wms?service=WMS&version=1.1.0&request=GetMap&layers=sde:prism_spring_index&styles=&bbox=-125.0,24.2,-66.9,49.9166666666687&width=768&height=339&srs=EPSG:4269&format=image%2Fpng&transparent=true

The problem is I added a second raster to the table for a different plant. How can I set things up so that I can pass a param in the url to geoserver to specify a "plant=lilac" parameter so that the appropriate raster is displayed?

I tried using (CQL_FILTER) like this but it doesn't seem to work

http://<myserver>/geoserver/sde/wms?service=WMS&version=1.1.0&request=GetMap&layers=sde:prism_spring_index&CQL_FILTER=plant=%27lilac%27&styles=&bbox=-125.0,24.2,-66.9,49.9166666666687&width=768&height=339&srs=EPSG:4269&format=image%2Fpng&transparent=true

I'm pretty new to geoserver, so I'm not sure if I need to define two layers pointing to the same postgis table to get this working or if I just need one layer. Either way, how do the layers know how to specify which raster to pull from a table?

Here's a screenshot showing how CQL is disabled when I preview the layer. Even though you only see one US raster, there are multiple US rasters in the table that the layer points to and I'm not sure how it is picking which one to show. Is there a way to enable the CQL filter? I'm using imagemosaicJDBC.


texts: How to get most accurate Point-Point distance queries in Europe with Spatialite

texts:
I have a GeoDjango app, currently on PostGIS, which I am converting to use Spatialite db instead.  I have no specialised knowledge about GIS systems :(
The locations table has a Geography column storing a single point in SRID 4326 per row, originally derived from lat/long data.
It's a long time since I built the original app but I believe I had read advice like:
https://docs.djangoproject.com/en/3.1/ref/contrib/gis/model-api/#geography-type
https://postgis.net/docs/using_postgis_dbmanagement.html#PostGIS_GeographyVSGeometry
...and this led me to believe that I would get more accurate point-to-point distance measurements with out-of-the-box GeoDjango queries by using the Geography type.
I have loaded the data in Spatialite, the same table now has a Point column. The data was dumped and loaded via GeoDjango and the values look the same.
When I do the distance queries I get similar but different results between PostGIS and Spatialite. The discrepancy gets larger with larger distances.
If I inspect the method Django is using in PostGIS:
SELECT ...
  ST_Distance("locations"."geopoint", ST_GeogFromWKB(\'\\x0101000020e6100000c3bb5cc47762c0bf1893fe5e0ac14940\'::bytea))
    AS "distance"
  FROM locations
  WHERE ...

vs in Spatialite:
SELECT ...
  COALESCE(Distance("locations"."geopoint", GeomFromText(\'POINT (-0.128005 51.508129)\',4326), 0), 0)
    AS "distance"
  FROM locations
  WHERE ...

...this seems to fit with what I read about "Geography" vs "Geometry".
Would it be correct to say that the PostGIS Geography query is more accurate?
And reading here https://www.gaia-gis.it/spatialite-3.0.0-BETA/spatialite-cookbook/html/srid.html it sounds like the key to getting accurate distances in Geometry queries is to use an appropriate projection, would that be correct?
If so, is there a way to get more accurate distance queries in Spatialite by re-projecting my points in a different SRID?
My use case would be "Select all rows where distance from <some point> is less than X kilometres".  All the points are in Europe.


texts:
Has anyone used / come across MapMaker GIS Software http://www.mapmaker.com?

A consultant has stated that they will use this for a piece of work involving GIS analysis and that it is compatible with ArcGIS.

The work will involve analysing many datasets, predominantly in SHP format to produce a land management strategy - the outputs of which will need to feedback into our corporate system (PostGIS server / ArcMap desktop).

From the website it seems that it will import SHP files, but no mention of export. It looks like a cartographic tool rather than capable of GIS analysis.

I would be grateful if anyone has used it and can comment on suitability for task and interoperability.

Thanks


texts:
I've written an SQL query which works in the query section of the QGIS DB Manager.

SELECT Max(S_order), SegmentID, S_Order, Length, Buff_Segme, Buff_Order, Buff_Lengt, Buff_Uniqu
FROM "Channel_Network_Within_Buffers_With_SSSI_Ref"
GROUP BY "Channel_Network_Within_Buffers_With_SSSI_Ref"."Buff_Uniqu"

Now, I want to use that query in the Execute SQL tool as part of a process model. I've copied the exact same SQL query into Execute SQL but I get the following error:

  Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS3~1.6/apps/qgis/./python/plugins\processing\algs\qgis\ExecuteSQL.py", line 160, in processAlgorithm
  vLayer.fields(), vLayer.wkbType() if geometry_type != 1 else 1, vLayer.crs())
  Exception: unknown
  Execution failed after 0.28 seconds

This is how the Execute SQL query is set up:

Why would the same SQL query work in the DB Manager query builder but fail in the Execute SQL Tool?


texts: ArcPy - What's the syntax for comparing a variable against two values for Boolean?

texts:
I'm a student trying to get a Boolean True/False value based on if a raster cell value falls within a given range, and I want to clean it up so that I can process each separate range as only one new raster. Basically, I just want to get values that are, say, between 5 and 20 for one raster, the values that are between 150 and 270 for another raster, and between 41 and 43 for a third raster, and then get the composite of those three.

The way I've got it written is fine, but I want to clean it up a little bit. However, I have two Boolean operators for each individual range, and every time I try to combine the operators into one per raster, ArcGIS tells me that it cannot determine the Boolean value, and that my statement is ambiguous.

Why does this happen? Here's what I have:

from arcpy.sa import *

outraster1 = arcpy.sa.Slope("Elevation")
outraster1.save("slope")
slope = arcpy.Raster("slope")

outraster2 = arcpy.sa.Aspect("Elevation")
outraster2.save("aspect")
aspect = arcpy.Raster("aspect")

landuse = arcpy.Raster("landcover.tif")

goodslope = 5 < slope
goodslope2 = slope < 20
bestslope = goodslope & goodslope2

goodaspect = 150 < aspect
goodaspect2 = aspect < 270
bestaspect = goodaspect & goodaspect2

goodforest = 41 <= landuse
goodforest2 = landuse <= 43
bestforest = goodforest & goodforest2

bestfinal = bestslope & bestaspect & bestforest
bestfinal.save("challenge9_1")

The above gives me what I want, but I'm not sure why I can't just use this:

goodslope = ((slope>5) and (slope<20))
goodaspect = ((aspect>150) and (aspect<270))
goodforest = ((41 <= landuse) and (landuse <= 43))
goodfinal = goodslope and goodaspect and goodforest
goodfinal.save("challenge9_1")

Is there a reason why the syntax above keeps giving me an ambiguous error?


texts:
I have a WGS84 coordinate pair (37.7749295, -122.4194155). I have an image that is a Robinson Map Projection of the world. I'm trying to write a python script to locate my WGS84 coordinate on my Robinson Map.

As an essential step in this, I'm trying to use pyproj, the python-adaption of Proj4, to convert my WGS84 coordinate to a Robinson coordinate. Unfortunately I am new to Proj4, and somewhat confused by the documentation. I hope that someone can help me out.


texts: How to get extent from current resolution in OpenLayers 3 for printing

texts:
I need to print a map with a fixed scale of 1:500 horizontally on a A4 paper (297x210mm). The map is done with Openlayers, the print request is sent to QGIS server using WMS getPrint function.

Till now I used a variable scale, printing the map in the actual visible resolution:

var dpi = 120;
var a4_width = 297;
var a4_height = 210;
var inch = 25.4;

var width = Math.round(a4_width * dpi / inch);
var height = Math.round(a4_height * dpi / inch);

var extent = map.getView().calculateExtent([width, height]);

window.open("http://localhost/cgi-bin/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetPrint&FORMAT=pdf&TRANSPARENT=true&LAYERS=all&CRS=EPSG:3857&map0:STYLES=&map0:extent="+extent+"&TEMPLATE=DIN_A4_horizontal (1:500)&DPI=120");

How can I calculate the extend of the printing area using the center part of the actual view?

Probably the solution is part of this answer, but I can't figure out how.


texts:
I have a set of points in 2 dimension plan : A(x1,y1), B(x2,y2), C(x3,y3) ....

How can I know which point of these points (A,B,C, ...) is the centroid of the other points


texts: How can I see the overlay of two features (point & polygon)?

texts:
I am trying to open two digitized layers, one containing points in a region and second which I have downloaded separately for the region as polygon layer. I would like to overlay two of these layers to see which point belongs to which area.
The problem is, I can not see an overlay of layers (points on polygon).


texts: Using Well-Known Text and geocoding on Load

texts:
So, I have tested adding columns with the WKT format into Carto, and it has always worked great (I tried 1 row, then 1000 random rows from our data). I geocoded it on load by naming the column with the WKT the_geom. When I uploaded all the data, some of the data was corrupted (random numbers and letters, no format). After testing, I have found the following:

The csv being used is normal and fine.

When the data is mixed (some the_geom values are blank, some are in WKT format), less than
40,000 rows can be added (40K doesn't work, 39.5K does)

The data that is corrupted is the first 24,000-25,000 rows (it varies) that have a value in the_geom column (only those first ones, no matter whether its 40K rows or 1.1 million).

If I upload 39K rows, the same rows that were corrupted before are now fine and geocoded.

I can load as many rows as I want, as long as the_geom is in WKT format OR blank, but not mixed.

All of this to say, I found a work around (upload 2 separate tables [one with the_geom WKT format, the other blank], then using sql union to combine them). Unfortunately, I cannot make a table out of these 2 because I have reached the size limits for our account (Xentity). Any ideas what is causing this? If not, is there any way we could do a one time extension of the DB size, in order to combine the tables (then delete the 2 halves)?


texts: Setting max/min values of Singleband grey layer using PyQGIS

texts:
I have a raster layer in QGIS and I would like to write Python code to change it's max and min values that would normally be under Layer Properties> Symbology> Band Rendering. The layer has a Singleband gray render. This is what I've written so far:
rLayer = iface.activeLayer()

rLayer.renderer().setClassificationMin(-426)
rLayer.renderer().setClassificationMax(8848)

Sadly it doesn't work for Singleband gray renders, only for Singleband pseudocolor.
Could anyone tell me what changes I should make to the code to make it work for the grayscale (Singleband gray)?


texts: Benchmarking QGIS and ArcGIS for Desktop speed for "significant buffer on very large shapefiles"?

texts:
I am a newcomer to the world of GIS. I realized that QGIS was relatively similar to ArcGIS for Desktop.

I would like to know if the two programs are comparable in terms of computational speed.

For example, will there be differences in computation time if I made a significant buffer on the same very large shapefile using QGIS and ArcGIS for Desktop?


texts: How to get features from a Service Layer?

texts:
I am using a Map Service from http://watersgeo.epa.gov/arcgis/rest/services/NHDPlus_NP21/WBD_NP21_Simplified/MapServer/. It offers multiple layers and its a great resource.

I'm wondering if its possible to access or create features layers from this service? In the end, I'd like to shade the watersheds according to differing data values similar to the example at https://developers.arcgis.com/javascript/jssamples/styling_svg_fb.html.

From what I'm reading at https://developers.arcgis.com/javascript/jshelp/inside_feature_layers.html, it seems like Service Layers just do not provide all the fun vector information that Feature Layers do. So, even if I can't get access to the color of the different shapes, can I at least get the gps coordinates so I can draw new polygons to work with?


texts: How can I add a new field in a raster attribute table and loop through it?

texts:
I am working with an existing raster attribute table.  I wonder how can I add a new field to it and then loop through the existing records to populate this new field? Of course I have done this many times for a feature class and a database table but I do not know how to do it with a raster table. I work with VBA. Any code please?
Thanks.


texts: ArcObjects - Create a map in the background from layers and save as JPG

texts:
Given a series of map layers is it possible to arrange them in some order, add basic map features (north arrow, scale bar) and then save the result as a JPG. I'd like to accomplish this in the background as a button-triggered event in an ArcMap Add-On.

I've seen ExportActiveView, but I think this will only work if I have the map layers already loaded into the ActiveView and on the users' screen.

Is it possible to create this map behind the scenes?


texts:
On GitHub is a wonderful installation instructions, but most users GitHub, use Linux, it's not worked on Windows....
https://github.com/openlayers/ol3/wiki/Developer-Guide

How can I build OpenLayers 3 on Windows?

Can't install Pystache

Traceback (most recent call last):
  File "C:\Openlayers\pystache-0.5.3\setup.py", line 111, in <module>
    import setuptools as dist
ImportError: No module named setuptools

And Closure Linter

Warning (from warnings module):

     File "C:\Openlayers\lib\distutils\dist.py", line 267
        warnings.warn(msg)
    UserWarning: Unknown distribution option: 'install_requires'


texts:
How to start learning javascript for GIS?

Should learn JS first and then try to use it for gis or there is any way to learn JS for GIS from beginning?

What is the site for learning js for GIS mapping?


texts: Create marker icon based on HTML code for OpenLayers

texts:
Is there any way that I can create a marker in OpenLayers with an HTML icon?

Example in Mapbox:

 marker.setIcon(
       L.divIcon({
            className: 'map_primary_marker_container',
            html: "
TEXT

",
            iconSize: [160, 50],
            iconAnchor: [80, 40]
       })
 );


texts: How to write a list of Polygons to a .shp file in Python

texts:
I create a list of Polygons and MultiPolygons outmulti as follows:
for pol in polygonToMulti:
    for pol2 in tre:
        if pol.intersects(pol2)==True:
            counter = counter+1
            # If they intersect, create a new polygon that is
            # essentially pol minus the intersection
            nonoverlap = (pol.symmetric_difference(pol2)).difference(pol2)
            outmulti.append(nonoverlap)

finalpol = MultiPolygon(outmulti)

I now want to save these Polygons and MultiPolygons of the outmulti list in a shp file & then be able to plot and visualize all the objects.
Firstly, I thought of keeping the values of wkt into a different list like this:
list = []
for i in range(len(outmulti)):
    list.append(outmulti[i].wkt)

and then I thought of creating a GeoDataFrame having as the geometry attribute the above list list.
gdf3=gpd.GeoDataFrame(geometry=list)
however this gives a TypeError saying it's not a valid geometry object.
Secondly, I tried creating a GeoDataFrame with a column name "geometry" and then writing all the wkt values of outmulti into that column and then saving the final GeoDataFrame into a file like:
#Make list of Polygons -> Geo Data Frame
gdf3=gpd.GeoDataFrame(columns={"geometry"})

for i in range(len(outmulti)):
    temp = gpd.GeoSeries(outmulti[i].wkt)
    gdf3= gdf3.append({'geometry':temp}, ignore_index=True)

    gpd.geometry

gdf3.to_file(filename='dPolygons.shp', driver='ESRI Shapefile')

But in that case I get the error:

AttributeError: 'Series' object has no attribute 'has_z'

I guess this is because the shp file must have some specific attributes in order to be created, but I cannot understand which and where to get those from.


texts: 'Selecting by freehand' tool - selection disappears

texts:
I have been trying to select some features of my map on QGIS to create the boundaries of a domain, however when I draw a shape using the 'Select by freehand' tool, the selection disappears as soon as I release the mouse.

Basically, I can only see the selection while I click and drag across the map, but I am unable to save what I have selected.

Can somebody help me figure out why?


texts: Is the native ArcGIS database capable of doing generalized data manipulation?

texts:
I'm doing political analysis across redistricting cycles, and need to do a lot of data manipulation. I need to split up precincts, do updates, joins, calculations, etc... I'll need a method that is able to handle a table of 500k census blocks. I just bought the home version of ArcGIS.

Should I try to store the data and do the manipulation using the internal .gdb file? That seems the easiest, but I have no idea if it's possible to do SQL-like updates and manipulations. I haven't been able to find the answer in the docs.

I'll need to do:

update table set this = that
create table as (select * from table1, table2 where someclause)

If ArcGIS is not capable of this, what do most people do instead?


texts: Satellite data between 1975-1985 in Brazil?

texts:
So I've been looking for satellite imagery to use for a project without any luck. I need images between 1975 and 1985 near the Roosevelt river in Brazil, from any satellite at any spatial/spectral resolution (could be aerial also!).

I've checked USGS Earth Explorer and there are only 4 images from Landsat 2 and Landsat 4, and they're all cloudy.

So I was just curious if anybody knew of a satellite/aerial imagery source that was flying around at that time.


texts: Altas with different numbers of page and maps for every feature

texts:
I want to create an atlas and the features to show are OUTSIDE the feature of the layer coverage.
The atlas feature 1 is Switzerland, but I want to show all embassy that are geographically OUTSIDE Switzerland, but that belongs to Switzerland, namely the Swiss embassy around the world. To show them I will create many maps and manually define the maps views. The same thing I will 2 with atlas feature 2 that is France and so on. To show them I need to create many different maps and pages for every atlas feature because every country (namely atlas feature) has different numbers of embassy and so on.
Switzeraland has x embassy around the world, French has y embassy and so on. So every atlas feature need different numbers of page or maps shape inside the page. My problem is that I don't know how to create an atlas with different number of pages and map shapes. It seems me that atlas can work only recurrent. If feature 1 has 3 page, feature 2 can have only 3 page, and feature 3 can have only 3 page. How can I have an atlas that change the page number for every atlas feature?
I add an example. You see an atlas with maps shape differents and in feature 1 you see the there are two page, in feature 2 there is 1 page. So, feature 1 has not the same layout of feature 2 (one vs two page and so on)

Someone could ask me, why do yo need an atlas? I need an atlas because it helps me to automatically create a filter of the vector embassy that otherwise would have shown the embassy of the entire world. Feature 1 of the atlas that is Switzerland allows me with the filter in the symbology to see automatically only the Swiss embassy from the vector embassy and feature 2 only from France and so on. Atlas has only this purpose.
Pay attention: I don't work with embassy but I used only as example to show the concept in this thread.


texts: ArcGIS Javascript API Editing - Create Feature from Geometry Object

texts:
I'm interested in passing a geometry object with certain attributes to an exposed feature class for committed edits. Essentially, I want to manually create a new feature without the template pane or editor widget. The geometry and attributes are both pre-defined, and the user just needs to select one pair from a group and push an "Add to Geodatabase" button.

I'm reading through the documentation and I can't see a clear way to do this. I see that if I don't want the pane or widget, then I need to use a geometry service, but I don't see a "create feature" method, or something that accepts a geometry object or attributes.


texts: With QGIS, why use QT Designer tools for creating Slots/Signals/ etc, if plugin coding can do it alone

texts:
I am curious why there is facility within QT Designer to create signals and slots as, to my surprise, I did not need to do this.   All of the necessary coding was accomplished completely within the plugin, and without using QT Designer.

I installed QGIS 2.6.0, lately, and have successfully created a signal / slot in a test plugin (an effort beginning with PluginBuilder): getting a ListView widget to respond to a combobox change

Can someone explain to me a) what is the use of the slot/signal facility in QT Desginer, when using QGIS and python and b) how is that slot/signal facility in QT Designer used?

Maybe a pointer to a full web-explanation on this point would be helpful.   The QT tutorials I've found, so far, don't help, or, at least, don't help in the context of a QGIS plugin.


texts: Doors disappear from AutoCAD file when imported into QGIS

texts: I have a .dwg CAD drawing of a hospital floor. When I convert the file to .dfx and drag it into QGIS the doors are no longer visible. The walls show up fine.


texts: Display labels set from within QgsProcessingAlgorithm in QGIS 3.x?

texts:
I'm writing a simple QgsProcessingAlgorithm and I need to set some labels (to display the azimuths of the sides of a polygon).
I can set the labeling expression alright (the expression is properly set and after the QgsProcessingAlgorithm has run, if I open the properties dialog and simply click OK, the correct labels will be shown) but can't get it to display automatically.
I've commented some lines that I thought could do the trick but don't seem to do anything actually.

def processAlgorithm(self, parameters, context, feedback):
    source = self.parameterAsSource(parameters, self.INPUT, context)

    outputs = {}
    results = {}

    #Polygons to Lines
    alg_params = {
        'INPUT': parameters['INPUT'],
        'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
    }
    outputs['PolygonsToLines'] = processing.run('native:polygonstolines', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

    #Explode Lines
    alg_params = {
        'INPUT': outputs['PolygonsToLines']['OUTPUT'],
        'OUTPUT': parameters['OUTPUT']
    }
    outputs['ExplodeLines'] = processing.run('native:explodelines', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

    results['OUTPUT'] = outputs['ExplodeLines']['OUTPUT']

    layer = QgsProcessingUtils.mapLayerFromString(results['OUTPUT'], context)

    pls = QgsPalLayerSettings()
    pls.isExpression = True
    pls.fieldName = 'CONCAT(format_number(degrees(azimuth(start_point($geometry), end_point($geometry))) % 180, 2), \'°\')'
    #pls.enabled = True
    layer.setLabeling(QgsVectorLayerSimpleLabeling(pls))
    layer.setLabelsEnabled(True)
    #layer.commitChanges()
    #layer.triggerRepaint()

    return results

EDIT:
Although I'm primarily interested in solving the problem for the QgsProcessingAlgorithm, it is worth noting that after running the two algorithms manually, executing the following code (after selecting the Exploded layer) in the python console doesn't work either.

layer = qgis.utils.iface.activeLayer()

pls = QgsPalLayerSettings()
pls.isExpression = True
pls.fieldName = 'CONCAT(format_number(degrees(azimuth(start_point($geometry), end_point($geometry))) % 180, 2), \'°\')'
pls.enabled = True

layer.setLabelsEnabled(True)
layer.setLabeling(QgsVectorLayerSimpleLabeling(pls))
#layer.labeling()

#layer.commitChanges()
#layer.triggerRepaint()


texts: Reproject an aerial layer in QGIS and then open with ArcEngine?

texts: I have an Aerial layer (.SID file) which I would like to re-project in QGIS, then I would like to open it in ArcEngine.  I've figured out how to do this with vector layers in GGIS, its just a right click, "Save As", but my Aerial layer does not have a "Save As" option.


texts: One of the 3 clipped Landsat images does not overlap with 2 others in ArcMap

texts:
I am analyzing forest cover change using 3 time-series Landsat data. The whole scene (185 x 185 km) images overlap in ArcMap display completely. However, when I subset/clip the images to my study site using a polygon, one of three images does not overlap with other two. It shows a bit displacement (nearly half pixel). Hence I can not go further with my analysis of forest change change.

To note, all layers/images have same datum/projection. I tried to image-to-image georegister the non-overlapping image with another image to nearly 0 rmse accuracy. However, when I clipped/extracted the study area portion again from this image with that polygon, the image again was not overlapping with the other two!

How can I get all these images spatially overlap completely?

What are the possible problems which might arise in finding forest cover and biomass change for 3 time-series years using those 3 images?


texts:
I have drivetest data, and usually I do an interpolation using the values taken in the drive with their respective lat/long. I use natural neighbor. You can see an example attached.

I am trying to automate the process; however, in MapBasic, I cannot find any resource of a code to do that. I found something interesting which is vertical mapper, but it is like the GIS -> Grid Analysis in PlanetEV and I can't find a way.

Has someone tried to do this before?


texts: SQL check if first value in line is number

texts:
In the layer definition query in ArcGIS 10.4 I want to only use objects that look like

1/12 Main Street
12-2 Side Street
etc

and not ones that have

Main Street
Side Street

Is this possible?

So something like

Select *.FROM DCDB WHERE

    Address.startswith[0..9]


texts: Using .map() to display multiple images in an ImageCollection

texts:
How do I display multiple images from an image collection using .map()?

// Create clipping function
var clipper = function(image){
  return image.clip(geometry);
};

// Load image collection
var modis = ee.ImageCollection('MODIS/006/MOD09GQ').map(clipper);
var months = ee.List.sequence(1,12);

var byMonth = ee.ImageCollection.fromImages(
      months.map(function(m) {
        return modis.filter(ee.Filter.calendarRange(m, m, 'month'))
                    .filterBounds(geometry)
                    .median()
                    .set('month', m);
}));

// Create a visualization color scheme
var medVis = {
  min: 0,
  max: 8000,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};

var display = function(image){
  return Map.addLayer(image,medVis);
};

byMonth.map(display);

I get the following errors when I run the last few lines of code (please use any geometry of your choice): Please reopen the question as a user has a more relevant answer.

Can't encode an undefined value.

Layer 2: Layer error: Failed to decode JSON.
Error: Field 'value' of object '{"type":"ArgumentRef","value":null}' is missing or null.
Object: {"type":"ArgumentRef","value":null}.

Layer 3: Layer error: Failed to decode JSON.
Error: Field 'value' of object '{"type":"ArgumentRef","value":null}' is missing or null.
Object: {"type":"ArgumentRef","value":null}.


texts:
I would like to convert a GeoTIFF to PNG or JPEG RGB, and my GeoTIFF is stripped of a lot of valuable data (I think..)

I am on ubuntu and I am using gdal_translate to try and make a GeoTIFF turn into a regular PNG.

The GeoTIFF I have 'works'. I know because I can view it in QGIS on my Mac.

My GeoTIFF has 8 bands and I dont know which band corresponds to RGBA or what the other 4 bands stand for. It appears none of the bands have a valid 'color table'

  gdal_translate -of PNG test.tiff output.png

gives me

  ERROR 6: PNG driver doesn't support 8 bands.  Must be 1 (grey), 2 (grey+alpha), 3 (rgb) or 4 (rgba) bands.

When I try and just use any combination of bands for RGBA like this:

 gdal_translate -of PNG -b 1 -b 2 -b 3 test.tiff output.png

I just get a black image

In QGIS if I right click on my image and go to properties > Symbology, I see that the redband is band 1 (grey) (min 0, max 439, the green band is band 2 (min 0, max 460), the blue band is band 3 (min 0, max 454). and the image appears nicely in rgb - so the GeoTIFF 'works'.

I've also tried to make a gray image (even though I want an RGB), just because it looks like QGIS is saying band 1 is gray

  gdal_translate -of PNG -expand gray -b 1 test.tiff output.png

but that gives me

  ERROR 1: Error : band 1 has no color table

I dont really know what else to try, or if my image is black because I need to do something with 'scales' or something?

When I run tiffdump on my file I get

Magic: 0x4949 <little-endian> Version: 0x2a <ClassicTIFF>
Directory 0: offset 8 (0x8) next 0 (0)
ImageWidth (256) SHORT (3) 1<2604>
ImageLength (257) SHORT (3) 1<2233>
BitsPerSample (258) SHORT (3) 8<16 16 16 16 16 16 16 16>
Compression (259) SHORT (3) 1<1>
Photometric (262) SHORT (3) 1<1>
StripOffsets (273) LONG (4) 2233<18288 59952 101616 143280 184944 226608 268272 309936 351600 393264 434928 476592 518256 559920 601584 643248 684912 726576 768240 809904 851568 893232 934896 976560 ...>
SamplesPerPixel (277) SHORT (3) 1<8>
RowsPerStrip (278) SHORT (3) 1<1>
StripByteCounts (279) LONG (4) 2233<41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 41664 ...>
PlanarConfig (284) SHORT (3) 1<1>
ExtraSamples (338) SHORT (3) 7<0 0 0 0 0 0 0>
SampleFormat (339) SHORT (3) 8<1 1 1 1 1 1 1 1>
33550 (0x830e) DOUBLE (12) 3<2.20356e-06 2.20356e-06 0>
33922 (0x8482) DOUBLE (12) 6<0 0 0 -82.0121 27.3184 0>
34735 (0x87af) SHORT (3) 32<1 1 0 7 1024 0 1 2 1025 0 1 1 2048 0 1 4326 2049 34737 7 0 2054 0 1 9102 ...>
34736 (0x87b0) DOUBLE (12) 2<298.257 6.37814e+06>
34737 (0x87b1) ASCII (2) 8<WGS 84|\0>

gdalinfo returns

Band 1 Block=2604x1 Type=UInt16, ColorInterp=Gray
Band 2 Block=2604x1 Type=UInt16, ColorInterp=Undefined
Band 3 Block=2604x1 Type=UInt16, ColorInterp=Undefined
Band 4 Block=2604x1 Type=UInt16, ColorInterp=Undefined
Band 5 Block=2604x1 Type=UInt16, ColorInterp=Undefined
Band 6 Block=2604x1 Type=UInt16, ColorInterp=Undefined
Band 7 Block=2604x1 Type=UInt16, ColorInterp=Undefined
Band 8 Block=2604x1 Type=UInt16, ColorInterp=Undefined


texts: Preventing labels from overlapping certain features using QGIS

texts:
I have two layers in my QGIS Project (QGIS 3.4.8). One in the foreground, which contains surfaces, and one in the background, which contains a point with labels, see image below:

The problem is that the labels from the second layer are always in the foreground:

I would like the second layer's label to be displayed in the background, just like the layer itself.
Is this possible in QGIS?
Desired target result (image edited with MS Paint):


texts:
i have 26 pictures of the globe and i want to bound them. How can i proceed to abtain a flatworld map? it looks circa like this but in 26 pieces.
map

Can somebody help me?


texts: CartoDB: Trigger InfoWindow from inside another InfoWindow

texts:
Is there a easy way (maybe even generic solution) to trigger opening the InfoWindow of Marker 1 from inside the InfoWindow of Marker 0?

The goal would be to «jump» from marker to marker right from inside the InfoWindows (à la «Next Location»-Links).


texts: What are the effects of changing the default Staging folder path in ArcGIS Server 10.1?

texts:
I want to change Staging folder path in ArcGIS Server 10.1.

The reason is that in my C drive very low space is available and my default staging folder was created at

  C:\Users\\AppData\Local\ESRI\Desktop10.1\StagingFolder

and I want to change the location on D drive.Going through web help also. So are there any effects on ArcGIS Server connection settings and Existing map services?

Thanks in advance : )


texts: Using AddDataDialog in standalone ArcObjects/C# application?

texts:
I need some help with showing an "AddDataDialog" in a C# standalone application to select some feature classes from inside a file-geodatabase. Whatever I try, it ends up with an Exception which makes no sence to me. I used a ESRI snippet to provide the hWnd of  "this" application and a referenced bool.

Is it possible to use this Dialog out of the context of ArcMap? And, if yes, has someone out there an idea what I might have done wrong so far?


texts: Automatically cutting polygons created from Prepare Time Lapse Polygons tool to prevent polygon crossing a Network Restriction?

texts:
I am trying to develop isochrones (polygons showing how far you can get from a bus stop in a certain amount of time) for an entire transit network. In order to do this, I have created a custom network dataset that incorporates the transit system's GTFS and the pedestrian network using the Add GTFS to network dataset toolbox.

I have managed to develop accurate travel distances using my custom network datasets and the toolbox, but now I am hitting a snag with the automatically developed polygons.

Simply solving a Service Area in Network Analyst doesn't capture how far you can get at ANY time, so to get around this, I used the Prepare Time Lapse Polygons tool to show how far you can get in 15 minutes by transit at any time between 8 and 8:30am then dissolved all polygons with the same Facility ID. Unfortunately, the polygons do not account for barriers that would prevent a person from actually reaching a location in 15 minutes, ie. Highways.

The image shows an example of this issue. The green point is one bus stop that is located very close to a highway, shown as the red line. The pink polygon is the area that was developed by the time lapse polygons tool. The area outlined by the black dashed line is not actually accessible within 15 minutes because of the physical barrier of the highway. However, the polygon is showing it as accessible within 15 minutes.

I have attempted incorporating a network restriction in the Service Area prior to running the Time Lapse polygons tool, but this only prevents the lines generated from the Service Area from crossing the highway barrier.

My question is, is there a way to automatically prevent the polygons from showing these inaccessible areas as accessible, eg. cut out the area outlined in black dashes, and for any similar locations? Otherwise I will have to cut several hundred polygons manually.

example of polygon issue


texts: Representing GIS data in ArcScene so it looks like bar and circle?

texts:
Is there a way to represent GIS data in ArcScene  so it looks like a bar and circle?

I can only either get the bar by extruding my points but I wanted a circle on top. I don't know how to get my circles to float to the top of my bars.


texts: How to search on values of a value list in python toolbox?

texts:
I developed a simple python toolbox that lists values of a table in the value table parameter. The user can rank each row using the Rank field drop down values.Everything is ok untill user selected value. Now i want to update the table based on the rank values.My solution is search rank values in the value table (params3) and update params4 from these values.The main problem in search cursor is access to the Rank field ( search cursor is in execute function of my code). I think if i can search these values then i can use updatecursor to update the values in the table.

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Tool]

class Tool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""

        self.label = ""
        self.description = ""
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
         # line station parameter
        params0 = arcpy.Parameter(
        displayName="Line Station",
        name="line_station",
        datatype="GPFeatureLayer",
        parameterType="Required",
        direction="Input")
        params0.filter.list = ["Polyline"]

        # station parameter
        params1 = arcpy.Parameter(
        displayName="Station",
        name="point_station",
        datatype="GPFeatureLayer",
        parameterType="Required",
        direction="Input")
        params1.filter.list = ["POINT"]

        # Table
        params2 = arcpy.Parameter(
        displayName="Table",
        name="table",
        datatype="GPTableView",
        parameterType="Required",
        direction="Input")
        #  value table
        params3 = arcpy.Parameter(
            displayName='Values',
            name='values',
            datatype='GPValueTable',
            parameterType='Required',
            direction='Input')
         # Table output
        params4 = arcpy.Parameter(
        displayName="Tableoutput",
        name="tableoutput",
        datatype="GPTableView",
        parameterType="Derived",
        direction="Output")
        params4.parameterDependencies = [params2.name]

        params3.columns = [['Long', 'Station Code'],['GPString','Station Name'],['GPString','Rank']]
        params3.filters[2].type="ValueList"

        params = [params0,params1,params2,params3,params4]
        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    def updateParameters(self, parameters):
        if parameters[2].altered:
             if not parameters[3].altered:
                result = arcpy.GetCount_management(parameters[2].value)
                count = int(result.getOutput(0))
                with arcpy.da.SearchCursor(parameters[2].value,["Station","Line","Row"]) as cur:
                   vtab = []

                   for row in cur:
                        vtab.append([row[0],row[1],row[2]])
                        parameters[3].value = vtab

                parameters[3].filters[2].list = range(1,count+1)

        return

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

    def execute(self, parameters, messages):
        """The source code of the tool."""

        if  parameters[3].altered:

            with arcpy.da.SearchCursor(parameters[3].valueAsText,parameters[3].columns[2][1]) as curs:
                    vtab = []
                    for rows in curs:
                        vtab.append(rows[0])
                        # The code rised an error and can not find the Rank field.

        return


texts: How to Extract Thalweg line from DEM using ArcMap/ArcGIS Pro

texts:
I would like to use a DEM to extract the thalweg line from a river. I have bathymetric data of the Saguenay River.

As a final product I need my thalweg line as polyline. I'm using both ArcGIS 10.6 and ArcGIS Pro.

I tried it using this workflow and as final product I'm getting a line, but I am not sure that it's a thalweg line.

Fill tool (input is my bathymetric DEM)
Flow Direction tool
Flow accumulation
Con
Raster to polyline


texts:
I hope somebody knows this error message, since I managed before the downloading with the same script (last year):

source("ModisDownload.R")
mrtpath <- '/home/usuario/MRT/MRT/bin/'

modisProducts()

Product to download is MOD11A2, no 31

x <- 31

H <- c(10)
V <- c(10)

Date <- c('2016.02.01','2016.03.02')

try (ModisDownload(x=x,h=H,v=V,dates=Date,mosaic=F, proj = T, MRTpath = mrtpath, bands_subset="1 1 0 0 0 0 0 0 0 0 0 0 0",proj_type = "UTM", utm_zone=18, datum ='WGS84', pixel_size = 500))

producing this error:

***Error in strsplit(unlist(lapply(strsplit(items, "href"), function(x) { :
  argumento de tipo no-carácter***

This error is part of the script Modisdownload.R, which I didn't changed.
Any ideas??


texts: Seeking workaround for wrong identified prj files

texts:
I recently had a problem in QGIS which is - i guess - related to gdal and the identification of WKT definitions in the prj files. The problem is that the CRS ETRS89 / UTM zone 32N in Shapefiles coming from ArcMap are automatically (without message in QGIS) identified as ETRS89 / ETRS-TM32 - which is wrong.

The definition of ETRS89 / ETRS-TM32 is:

+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

And ETRS89 / UTM zone 32N is:

+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

So it seems to be the same, also the OGC WKT definitions are - as far as i see it - the same. But the ESRI WKT is - of course - different:

ESRI:

PROJCS["ETRS_1989_UTM_Zone_32N",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",
SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

OGC WKT:

PROJCS["ETRS89 / UTM zone 32N",GEOGCS["ETRS89",
DATUM["European_Terrestrial_Reference_System_1989",
SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4258"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],
PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],PARAMETER["false_northing",0],
AUTHORITY["EPSG","25832"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

So, my questions:

Is there any possibility to adjust the wrong identification via gdal, so that QGIS can identify and set the CRS correctly?
Is there any difference between those CRS?
And if not - why do they have two different EPSG codes?


texts: Is it possible to pass multiple WMS layer to a processing script

texts: Is it possible to pass multiple WMS layers as parameters to a QGIS 3 processing script, similar to QgsProcessingParameterMultipleLayers but using WMS or WFS Layer? QgsProcessingParameterMultipleLayers does only accept Raster-/Vectorlayers depending on the defined SourceType in QgsProcessing.


texts: Problems displaying joined attribute data through ArcGIS Server

texts:
I have an ArcSDE feature class in an ArcMap MXD which I have 'joined' to a view into an external SQL database table, that is connected via an OLE DB connection.

In the MXD the join works perfectly; I can see the joined attributes when I open the attribute table of the original ArcSDE feature class.

However when I then publish the MXD to a service on the ArcGIS Server and add it to a new application it only displays the original SDE feature class attributes and not the full 'joined' attribute dataset.

Anyone seen this before or have any ideas?


texts: could not find function "merge_eurostat_geodata"

texts:
Just a couple of codes with the same error at the time of merge data:

FIRST CODE:

library(eurostat)
library(dplyr)
library(ggplot2)
library(stringr)
fer <- search_eurostat("fertility")
get_eurostat("demo_r_frate2", time_format = "raw", stringsAsFactors = FALSE)  %>%
dplyr::mutate(cat = cut_to_classes(values, n=7, decimals = 1)) %>%
merge_eurostat_geodata(data=.,geocolumn="geo",resolution = "20",  output_class = "df", all_regions = FALSE) %>%
ggplot(data=., aes(x=long,y=lat,group=group)) +
geom_polygon(aes(fill=cat),color="dim grey", size=.1) +
scale_fill_manual(values =  c("#d73027","#fc8d59","#fee08b","#ffffbf","#d9ef8b","#91cf60","#1a9850")) +
# scale_fill_continuous(trans = 'reverse', ) +
guides(fill = guide_legend(reverse=T, title = "Total fertility rate (%)")) +
labs(title="Total fertility rate, by NUTS-3 regions, 2014",
subtitle="Average number of live births per woman",
caption="(C) EuroGeographics for the administrative boundaries
Map produced in R with a help from Eurostat-package  <github.com/ropengov/eurostat/>") +
theme_light() + theme(legend.position=c(.8,.8), text=element_text(family="opensans")) +
coord_map(project="orthographic", xlim=c(-12,44), ylim=c(35,70))

  Error in merge_eurostat_geodata(data = ., geocolumn = "geo",
  resolution = "20",  :  couldn't find function "merge_eurostat_geodata"

SECOND CODE:

# Data from Eurostat
current_pop <- get_eurostat("demo_r_pjangrp3", time_format = "raw",   stringsAsFactors = FALSE) %>%
filter(time == 2015,
     nchar(geo) == 4,
     sex == "T",
     age == "TOTAL") %>%
select(-unit) %>%
rename(current = values,
     time_current = time)
projected_pop <- get_eurostat("proj_15npms", time_format = "raw",  stringsAsFactors = FALSE) %>%
filter(time == 2050,
     nchar(as.character(geo)) == 4,
     sex == "T",
     age == "TOTAL") %>%
rename(projected = values,
     time_projected = time)

 left_join(current_pop,projected_pop) %>%
 mutate(change = projected / current * 100 -100) %>%
 dplyr::mutate(cat = cut_to_classes(change, manual = TRUE, manual_breaks =  c(-50,-10,0,10,25,130))) %>%
 # dplyr::mutate(cat = cut_to_classes(change)) %>%
 # merge with geodata
 merge_eurostat_geodata(data=.,geocolumn="geo",resolution = "20",  output_class = "df", all_regions = FALSE) %>%
 # plot map
 ggplot(data=., aes(x=long,y=lat,group=group)) +
 geom_polygon(aes(fill=cat),color="dim grey", size=.1) +
 scale_fill_manual(values = c("dim grey","#d7191c","#fdae61","#ffffbf","#a6d96a","#1a9641")) +
 # scale_fill_continuous(trans = 'reverse', ) +
 guides(fill = guide_legend(reverse=T, title = "%")) +
 labs(title="Projected percentage change of the population, by NUTS2 regions, 2015-2050 (%)",
   subtitle="%",
   caption="(C) EuroGeographics for the administrative boundaries
   Map produced in R with a help from Eurostat-package <github.com/ropengov/eurostat/>") +
 theme_light() + theme(legend.position=c(.8,.8),   text=element_text(family="opensans")) +
 coord_map(project="orthographic", xlim=c(-12,44), ylim=c(35,70))

And same error:

  Error in merge_eurostat_geodata(data = ., geocolumn = "geo",
  resolution = "20",  :  couldn't find function "merge_eurostat_geodata"

Any idea on how to solve it?


texts: Filtering polygons by long/lat values using sf package in R?

texts:
Suppose t is:

t <- structure(list(structure(list(structure(c(-89.990791, -89.990772,
-89.990901, -89.99092, -89.990791, 30.727025, 30.727083, 30.727114,
30.727057, 30.727025), .Dim = c(5L, 2L))), class = c("XY", "POLYGON",
"sfg")), structure(list(structure(c(-89.991691, -89.991755, -89.991755,
-89.991691, -89.991691, 30.716004, 30.716004, 30.715916, 30.715916,
30.716004), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"
))), class = c("sfc_POLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = -89.991755,
ymin = 30.715916, xmax = -89.990772, ymax = 30.727114), class = "bbox"), crs = structure(list(
    epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L)
> t
Geometry set for 2 features
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -89.99176 ymin: 30.71592 xmax: -89.99077 ymax: 30.72711
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POLYGON ((-89.99079 30.72703, -89.99077 30.7270...
POLYGON ((-89.99169 30.716, -89.99175 30.716, -...

How can I filter polygons by long/lat boundaries? Suppose I want to filter out any polygon that hax lat>30.72 (so to keep only the second polygon). Is there any specific function I can use to filter polygons?


texts: Does OSM2PO take into consideration turn restrictions?

texts: It is supposed that for using Shooting star function in pgrouting, we need to have RULE and TO_COST fields in our table, but OSM2PO hasn't created those fields. what must i do for OSM2PO takes into consideration the turn restrictions?


texts: Does GeoServer have Public directories reachable over internet?

texts:
I am trying to figure out whether or not there are any directories within GeoServer latest which allow me to offer image files over the internet. I am using windows server 2012. I have internet access to the endpoints (set up through a reverse proxy) which allow users to see the map services e.g servername/arcgis/services/directory/servicename/MapServer/ but as far as I know these are virtual directories and do not allow me to save images within them.

I am looking to simply find an open directory without having to set up another reverse proxy.


texts: Is it a bad idea to hold references to ArcObjects inside business classes?

texts:
According to best practices is it a bad idea to cache features and geometries inside business objects? If it is bad, what is the alternative way?

References to features may be replaced with their Ids, but what about geometries?


texts: Inserting dynamic text: way to set text format and placement preferences?

texts:
I have a series of maps that I need to insert a footer with the filepath and mxd name.  I see that I can insert these using dynamic text, which is great.  However, the text is always placed in the middle of the layout in a text format I do not like.

Is there a way to set up preferences so when you select "Insert\Dynamic Text\Document Path" the filepath is always placed with preferred settings(lower left, 6 pt font, etc)?

I am building the preferences into our templates so we will have this going forward.


texts:
I have downloaded a shapefile (.shp) that I'd like to read in using fiona.open and delete records. The reason I'd like to delete records is that some records have "None" geometries.

So I want to delete a record if it doesn't have a geometry: {'geometry': None, 'type': 'Feature', 'id':  etc etc etc}

Here's my code:

import fiona
with fiona.open('data/ssc_2016_aust_shape/SSC_2016_AUST.shp') as file:
    for p in file: # for every record - note that p is a dictionary
        if p["geometry"] is None: # if p doesn't have a geometry
            del p # then delete the entire record in the shape file

I'd like for the original shapefile (.shp) to be overwritten (this is important).

You can download the shapefile here: http://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/1270.0.55.003July%202016?OpenDocument

The file on the ABS website is called "State Suburbs ASGS Ed 2016 Digital Boundaries in ESRI Shapefile Format" in the form of a ZIP file.


texts: Batch maps (atlas plugin style) in same location

texts:
I've got more than 200 shapefiles of points representing different information for the same location. I would like to make as many maps as shapefiles with identical background and extent. I'd rather use QGIS, but ArcGIS 9.3/10 would be ok too.
I found EasyPrint and Atlas, which are great, but both are zooming on the shapefile.. which I don't want to. Atlas is way better as it lets you use your own print composer.. but it only wants one shapefile (coverage) while I've got many.

Any idea on what I could use? I don't mind a bit of programming, especially because that would save me a looooot of time!


texts: Adjust size of clustering icon in leaflet

texts:
In my leaflet application, I would like to reduce the size of cluster icons by half and also modify the font size of the value inside of the icon:

Having read the documentation on github, I have tried changing values under iconSize in the iconCreateFunction:

return new L.DivIcon({
  html: '
' + childCount + '
',
  className: 'marker-cluster' + c,
  iconSize: new L.Point(40, 40)
});

But this does not appear to change anything.  Is there a known workaround for changing the size of cluster icons and their text?


texts: Reduce region applied to multi-country area when ImageCollection is organized by country and year in Google Earth Engine

texts:
I am calculating zonal statistics for 0.01-degree tiles that fall within buffer areas across sub-Saharan Africa. I want to determine the estimated population counts within the tiles for years 2001-2012 using the Worldpop dataset. Earth Engine stores this image collection so that each feature (image) represents a country-year. When I am deriving Worldpop zonal stats. for just one country, I have used the following function before (in this example, for tiles ("grid") across Syria):
var pop = ee.ImageCollection("WorldPop/GP/100m/pop")
  .filterDate('2000', '2012')
  .filterMetadata('country', 'equals', 'SYR');

function newCollectionToImage(collection){
var stack = ee.Image(collection.iterate(function(img, prev) {
return ee.Image(prev).addBands(img);
}, ee.Image(1)));

stack = stack.select(ee.List.sequence(1, stack.bandNames().size().subtract(1)));
return stack;
}

var tch = pop.map(function(img){
var y = ee.String(ee.Date(img.get('system:time_start')).get("year"));
return img.set("year",y).rename(y);
});

var tchImg = newCollectionToImage(tch);

var pop_Syria =
tchImg.select(['2000', '2001', '2002', '2003', '2004', '2005', '2006',
'2007', '2008', '2009', '2010', '2011'],
['pop_00', 'pop_01', 'pop_02', 'pop_03', 'pop_04', 'pop_05', 'pop_06',
'pop_07', 'pop_08', 'pop_09', 'pop_10', 'pop_11']);

var Syr_gridpop = pop_Syria.reduceRegions(grid,ee.Reducer.sum());

But what if instead, I wanted zonal statistics from polygons in multiple countries? When I tried removing the metadata filter (that designates country), the resulting stack of images that this function produces is all zeroes. Here is my code.
Is there a way in which I can combine the country-level features in the Worldpop data by year (resulting in a global image) and then perform the zonal statistics?


texts:
I need to create a layer that contains all values between 6 and 9. <= 6 >9. All other values should be 0.
How can I do this with using the raster calculator?


texts: Transforming GCS_North_American_1983 to GCS_WGS_1984 using ArcGIS Desktop?

texts:
I posted Is there projection issue with Natural Resources Canada (NRC) Canada Atlas hydrology data? describing a shift in some of the NRC data (not all).

Apparently, this shift has to do with conversion between Geographic coordinates, latitude and longitude and Transverse Mercator Grid coordinates.

Can someone please provide a guide on how to do this transformation in ArcGIS Desktop?

I would like to transform that data to either Clarke_1866_Lambert_Conformal_Conic or
GCS_WGS_1984.

Source Data Coordinate System:

GCS_WGS_1984
Angular Unit: Degree (0.017453292519943299)
Prime Meridian: Greenwich (0.000000000000000000)
Datum: D_WGS_1984
Spheroid: WGS_1984
Semimajor Axis: 6378137.000000000000000000
Semiminor Axis: 6356752.314245179300000000
Inverse Flattening: 298.257223563000030000


texts:
I'm wondering if this is the right way to do this, or if there's an easier way.

For illustration sake, take this SpatialPolygons object:

main = cbind(c(0, 0, 1, 1),
             c(0, 1, 1, 0))
hole = main/3 + 1/3
island = cbind(c(1.05, 1.05, 1.55, 1.55),
               c(0, .5, .5, 0))
P = Polygons(list(Polygon(main),
                  Polygon(hole, hole = TRUE),
                  Polygon(island)), ID = 1)
SP = SpatialPolygons(list(P))
plot(SP, col = 'red')

I'd like to clip the "island" polygon while retaining the topology of the "main" polygon -- namely, that is has a hole (in my case, it has several holes, but this should do).

Here's what I came up with:

#extract sub-polygons of this polygon
subP = SP@polygons[[1L]]@Polygons
#use the area to identify the largest
areas = sapply(subP, slot, 'area')
#keep holes or the largest actual polygon
keep = sapply(subP, function(p) p@hole | p@area == max(areas))
#construct a new SpatialPolygons object from scratch using index
rem = SpatialPolygons(list(Polygons(subP[keep], ID = 1)))
plot(rem, col = 'red')

That strikes me as being a bit circuitous -- am I missing a method that would handle this more directly? Other related questions I've seen speak more about a SpatialPolygons object that's made up of many polygons instead of one polygon and many Polygons (alas rendering this nuanced search all but impossible).


texts: Export dimensions of NetCDF as individual rasters in ArcMap

texts: I have a NetCDF file in ArcMap that I created using the Make NetCDF Raster Layer tool. The band dimension of interest is the variable "time", which has 316 grids. I need each individual grid as its own raster. However I cannot find any way to do this other than individually displaying each date, and exporting it to a raster. Is there a way to iterate through the grids and export each as a raster? Or am I stuck manually exporting each and every one?


texts: How to write it on ArcGIS Raster Calculator?

texts:
I just started using the Raster Calculator, so I don't know how to build the formulas.

How can I do this?

I need to write both of these formulas:

D = is a raster file called "Slope_2"

P = 13


texts: How to clip raster image by the extent of each row in a shapefile with ArcGIS 9.3

texts:
I'm using Python to clip a raster image by each row in a shapefile with ArcGIS 9.3.  I am hoping to end up with one raster image for each polygon row, clipped to the extent of the polygon row. However, my current script clips the raster to the entire shapefile, instead of setting the extent by each polygon row.  My script is as follows:

# Import system modules
import sys
import string
import os
import arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create(9.3)

# Load required toolboxes...
gp.AddToolbox("C:/Program Files (x86)/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")

# Designate a workspace
gp.workspace = "L:\Clip"

# Assign variables...
clipto = "test.shp"

# grids to clip to
clipit = "vfcm.tif" # features to be clipped

# Cursor time
rows = gp.SearchCursor(clipto)
row = rows.Next()

# clip management
while row:
    n = str(row.VAHU6) # assign a variable for the processing message based on a field
print "clipping to: " +n # tells what grid was clipped

# Clip_management <in_raster> <rectangle> <out_raster> {in_template_dataset}   {nodata_value} {NONE | ClippingGeometry}
outputname = str(row.VAHU6) + "_c"
desc = gp.Describe(clipto)
extent = str(desc.Extent)
gp.Clip_management(clipit, extent, outputname, clipto, "", ClippingGeometry) #update str if necessary
print str(row.VAHU6)
row=rows.next()

# reset the array...
del rows
del gp

Any ideas how to fix it?

Thanks!


texts: ArcGIS Javascript Print Text Symbol Color Change

texts:
I have an ArcGIS Javascript Map Application that uses the Print task to export the map to a PDF.  I have given the users the ability to add graphics to the map, such as points or polygons, along with the ability to add a TextSymbol graphic in order to customize the PDF map output.  The problem I am having is with the TextSymbol appearance on the PDF printout.  For some reason I cannot change the color of the text on the output.  Even if the text on the web map is red, the color of the text on the pdf is black.  I have applied code to change the color of the textLayers by performing an esriRequest.setRequestPreCallback to intercept the JSON going to the Print service.  However, the text still ends up black on the output.  My code for that is below.  I have done console.logs in order to see what the object looks like after changing the color and it does change it, but still the text ends up being black on the pdf output.  Does anyone know how to either change the color of the text on the PDF output or know how to add some type of rectangle background on the TextSymbol so that the labels can be more visible on the PDF outputs?

function changeTextColor(){
    esriRequest.setRequestPreCallback(function(ioArgs) {
        if (ioArgs.content && ioArgs.content.Web_Map_as_JSON){
            var webMapAsJson = ioArgs.content.Web_Map_as_JSON;
            var webMapObj = JSON.parse(webMapAsJson);
            for (var w in webMapObj.operationalLayers){
                if (webMapObj.operationalLayers[w].id = "map_graphics" && webMapObj.operationalLayers[w].url == undefined){
                    for (var l in webMapObj.operationalLayers[w].featureCollection.layers){
                        if (webMapObj.operationalLayers[w].featureCollection.layers[l].layerDefinition.name == "textLayer"){
                            for (var f in webMapObj.operationalLayers[w].featureCollection.layers[l].featureSet.features){
                                //currently I am just trying to change the text to be turquoise but it keeps returning black text on pdf output
                                webMapObj.operationalLayers[w].featureCollection.layers[l].featureSet.features[f].symbol.color[0] = 0;
                                webMapObj.operationalLayers[w].featureCollection.layers[l].featureSet.features[f].symbol.color[1] = 255;
                                webMapObj.operationalLayers[w].featureCollection.layers[l].featureSet.features[f].symbol.color[2] = 255;
                            }
                        }
                    }
                }
            }
            ioArgs.content.Web_Map_as_JSON = JSON.stringify(webMapObj);
        }
        return ioArgs;
    });
}


texts: Automate removal of duplicate polygon nodes from shapefile in ArcMap

texts:
I have a whole lot of polygon shapefiles, many of which have got duplicate nodes. This makes processing difficult as almost all tools (erase, union, calculate area) do not work.

I found that by removing duplicate nodes from the 'offending' polygons the shapefiles become usable. However there are well over a hundred shapefiles and going through all of these is a nightmare.

Is it possible to automate the process using ArcPy or VB? If it is possible, how can it be done?

If not, are there any free libraries in .NET (C# or VB.NET) or Python that can be used to perform this task?

I tried the Python shapefile library (http://code.google.com/p/pyshp) but I found that when it is done removing the duplicate nodes, ArcMap is still unable to handle the shapefiles, giving the same errors as before.


texts:
There is plenty unknown about how to use raster calculator, for example:

how to write NODATA (Null) value in expression
is it possible to use within (polygonLayerName) to clip result
between shown operators there is even "not equal" symbol, you need to find yourself it is "!=" etc.

Generally is there any manual/help/tutorial to QGIS raster calculator describing what is and is not possible to do with it?


texts: Dynamic postgis raster data type serving from geoserver

texts:
I have PostGIS database, which I am using to store raster image as raster datatype. A table has column date, raster, and other property. There are more than 30 tables in the database. The database tables are updated, where the raster for the particular time is replaced with the latest available and/or new rows are added.

I am doing little research on if I can use GeoServer to store the dynamic PostGIS DB and render as WMS/WFS for the web client. I need to be able to perform clipping, filter particular row based on date and style those data. What is the best approach to do this in geoserver? I followed postgis raster in geoserver, but it seems I need to generate pyramid and store them in a separate table using Image Mosaic JDBC. Since data are dynamic, is there another way without exporting the raster?


texts: pgrouting 2 pgr_createTopology not populating source and target field

texts:
I trying to integrate pgrouting 2.x. I am now using pgrouting 1.x. In version 1 I used assign_vertex_id. It already dropped in version 2 and became pgr_createTopology function

SELECT pgr_createTopology('roads_r', 0.00001, 'the_geom', 'gid')

It returning query result OK but its not populating source and target field of my road data. Is there any bug in this function?


texts: How to calculate path length of each point from the start using arcpy?

texts: As shown in the picture below I would like to calculate the path length of each point from the start. Can anyone give me some help in writing this in python for ArcGIS?


texts: Symbolize feature with most recent data if any exists

texts:
I am using ArcGIS Pro for this operation.  I have a attributes for features which represent a value for a given year.  A portion of the table is here:

As you can see, not every feature has a value for every attribute.  These values are the percentage of workers in any given year employed in an industrial sector job.  I would like to symbolize the features based on the most recent year.  That would mean that ArcGIS would use the 2011 value for Aruba, null for Afghanistan, null for Angola, the 2015 value for Albania, etc.

Is there a way to accomplish this task?


texts:
I have tons of LiDAR Bathymetric XYZ data of 5m resolution.
Data has a lot of noise/spikes which takes ages to clean for a small area. I have tried using Fledermaus for 3D cleaning the data. Is there any software available to automatically clean and visualize the data, using Standard Deviation?


texts:
I've just begun exploring QGIS relations option, now in version QGIS 2.2, details here: http://blog.vitu.ch/10112013-1201/qgis-relations

This is a great implementation of something I've been attempting to manually organize for a while now. Now, i'm wondering how to integrate the relations widget into a custom form (currently, it automatically adds only if using the auto generated form). I'm new to PyQt Designer, so i'm not sure where to start.


texts: QGIS - Expression to aggregate values from an array

texts:
In a calculated field in QGIS 3.10.2 I wan't to aggregate values from an array field. For this aggregation I'm using a filter to compare matching id's and return value from another field. I can aggregate values from a text or a numeric field, but I don't know how to reproduce this method with a list of values.
To be more precise:

myTargetLayer: a reference table which contains id and text values

idTargetLayer (integer)
myField

1
1
Grenache N

2
2
Cinsaut N

3
3
Carignan N

4
4
Viognier B

5
5
Chenin B

mySourceLayer: a layer which contains an array field of id's of target layer and a text field. In this text field I wan't to aggregate text values of the target layer (called 'myField' here). For example if array field is '1,2,3' I wan't to have 'Grenache N, Cinsault N, Carignan N'.

id
idTargetLayer (array)
Aggregate Field

1
9
4
Viognier B

2
8
1

3
7
1
Grenache N

4
6
2, 1

5
5
2, 1

6
4
1
Grenache N

7
5
1, 5

8
2
2, 1
Cinsaut N

9
19
2, 1
Cinsaut N,Cinsaut N

10
18
2, 1
,

11
17
2, 1
2,1

12
16
2, 1

Simple aggregation
aggregate(
layer:=layer_property('myTargetLayer','name'),
aggregate:='concatenate',
expression:="myField",
filter:= attribute(@parent,'id_source_layer') = attribute($currentfeature,'id_target_layer')
)

If I try to iterate through each value of the array, it returns the same duplicated value.
array_to_string(
    array_foreach(attribute($currentfeature,'id_source_layer'),
        aggregate(
        layer:=layer_property('myTargetLayer','name'),
        aggregate:='concatenate',
        expression:="myField",
        filter:= @element = attribute($currentfeature,'id_target_layer')
        )
    )
)

How to write an aggregate expression which iterates through each value of an array?


texts: Pixel shifts from its actual position in Matlab geospatial plots?


texts: Difference between SpatialJoin with "...where..." and "join...on..."

texts:
I wonder where the difference between these two spatial joins is and, if there is a difference, when to use which.

SELECT *
FROM points p, shapes s
WHERE ST_Within(p.geom, s.geom);

and

SELECT *
FROM points p
INNER JOIN shapes s ON ST_Within(p.geom, s.geom)

Is one faster in certain situations?


texts:
I tried to clip a DEM raster by extent, but I left it there for about 2 hours and got a 45GB file even though the original raster is actually 77MB. Since the process wasn't complete, I couldn't open the file. Based on what I have read, this might be related to the original raster being compressed... but could it really be this different (full data = 77mb, clipped region of interest = 45gb)? Also, I know that DEM files carry a lot of data, so it should take a decent amount of time... Should I leave it here for hours and forget about it?
I thought my computer could be the problem, but I have never had such a problem using QGIS.

QGIS Version: 2.16.3

Processor: 2.5GHz

Memory: 8GB RAM

Storage: 256GB SSD (with 100GB free)

This is the setup i tried.

Raster info:

  Driver: GTiff/GeoTIFF   Size is 5400, 3600   Coordinate System is:
  PROJCS["SIRGAS 2000 / UTM zone 23S",
      GEOGCS["SIRGAS 2000",
          DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",
              SPHEROID["GRS 1980",6378137,298.2572221010002,
                  AUTHORITY["EPSG","7019"]],
              TOWGS84[0,0,0,0,0,0,0],
              AUTHORITY["EPSG","6674"]],
          PRIMEM["Greenwich",0],
          UNIT["degree",0.0174532925199433],
          AUTHORITY["EPSG","4674"]],
      PROJECTION["Transverse_Mercator"],
      PARAMETER["latitude_of_origin",0],
      PARAMETER["central_meridian",-45],
      PARAMETER["scale_factor",0.9996],
      PARAMETER["false_easting",500000],
      PARAMETER["false_northing",10000000],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]],
      AUTHORITY["EPSG","31983"]]   Origin = (-48.000000000000000,-23.000000000000000)   Pixel Size =
  (0.000277778481481,-0.000277777777778)   Metadata:
  AREA_OR_POINT=Area     TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=IDRISI KILIMANJARO     TIFFTAG_XRESOLUTION=87.699997
  TIFFTAG_YRESOLUTION=87.699997   Image Structure Metadata:
  INTERLEAVE=BAND   Corner Coordinates:
  Upper Left  ( -48.0000000,
  -23.0000000) ( 45d 0' 0.00"W, 90d 0' 0.00"S)
  Lower Left  ( -48.0000000, -24.0000000) ( 45d 0' 0.00"W, 90d 0' 0.00"S)
  Upper Right ( -46.4999962, -23.0000000) ( 45d 0' 0.00"W, 90d 0' 0.00"S)
  Lower Right ( -46.4999962, -24.0000000) ( 45d 0' 0.00"W, 90d 0'
  0.00"S)
   Center      ( -47.2499981, -23.5000000) ( 45d 0' 0.00"W, 90d 0' 0.00"S)   Band 1 Block=5400x1 Type=Float32, ColorInterp=Gray
  Min=-0.000 Max=1.000      Minimum=-0.000, Maximum=1.000, Mean=0.639,
  StdDev=0.233     Metadata:
      STATISTICS_MAXIMUM=0.99998646974564
      STATISTICS_MEAN=0.63865509358276
      STATISTICS_MINIMUM=-0
      STATISTICS_STDDEV=0.23343316197389


texts:
I would like to select records 0-10 in INPUT. Algorithm works in QGIS Python console, but in cmd I'm getting error:
"Could not load source layer for INPUT:
  C:/Users/ja/Inne/Desktop/budynki/budynkishp.shp|subset=\"fid\" >=0 and \"fid\" <= 10 not found."

What is the problem?
processing.run("native:joinattributesbylocation", {'INPUT':r'C:/Users/ja/Inne/Desktop/budynki/budynkishp.shp|subset=\"fid\" >=0 and \"fid\" <= 10',
                                                   'JOIN':r'C:/Users/ja/Inne/Desktop/adresy/adresyzukladem.shp',
                                                   'PREDICATE':[1,5],
                                                   'JOIN_FIELDS':['PNA','SIMC_id','SIMC_nazwa','ULIC_nazwa','Numer'],
                                                   'METHOD':0,
                                                   'DISCARD_NONMATCHING':True,
                                                   'PREFIX':'',
                                                   'OUTPUT':r'C:/Users/ja/Inne/Desktop/caloscpython.shp'})


texts: Update Oracle DB table from Postgres/PostGIS DB query

texts: I have an Oracle DB where Geolocation (Longitude, Latitude) data is storing. I want to update another column which is the neighborhood of that geolocation. Currently, I do not have oracle spatial support. I have that neighborhood polygon stored in Postgresql DB in PostGIS format. I can easily update that geolocation to the neighborhood in Postgresql through ST_Contains function. Is there any way to update Oracle neighborhood column using PostgreSQL DB?


texts:
I am looking for ways to learn querying for SQL Server Spatial (2008/2012).  I have a table created from Shape2SQL and have been trying to run a simple intersection query to get results from that table that fall within a polygon I am creating.

I got as far as:

declare @poly geometry = 'POLYGON((-9372539.231 4862367.224, -9370110.628 4925510.907, -9457540.343 4926118.058, -945571889 4840509.795, -9372539.231 4862367.224))'

I haven't figured out how to take the geometry from my imported shapefile table and get the points that fall within the polygon geometry.  The options I have found online seem to take WKT from a field in the table to run the intersection.

I am looking for resources that apply to this kind of querying, or any other suggestions that may help me out along the way.  I picked up a copy of Beginning Spatial with SQL Server 2008 but haven't read through it yet.

Eventually, my goal is to be able to generate a query based off a map's extent and query features that fall within that extent.

Edit

I am looking for a query to allow me to do the above intersection.  Any additional resources would be great but not an answer for me.

So far, I have found how to grab the geom records and convert them to WKT.

SELECT [geom].STAsText() FROM [Table] WHERE [geom] IS NOT NULL

What I am trying to do is take that data set, currently I am putting it into a temporary table, but I don't think that is necessary.  What I want to do is find out which records from that data set fall within the generic polygon (@poly) I am setting up.

I am thinking if I take the WKT result and put it into a GEOMETRYCOLLECTION, I can pass that into something like SELECT @myWKTResults.STIntersection(@poly).ToString().  I just don't know/can't find a good example of the next step.

Result

Going from mapBaker's comments, I was able to set this up to work.  My query ended up being:

declare @poly geometry = 'POLYGON((-9486683.581 4810152.256, -9282073.762 4821688.121, -9262037.786 4625578.413, -9477576.319 4628614.167, -9486683.581 4810152.256))'

SELECT * FROM [TABLE]
WHERE @poly.STIntersects([geom].STAsText()) <> 0;

The only thing I really adjusted as adding in the .STAsText() into the query since running without didn't crash, but I also got 0 rows back.


texts:
I have a large set of 4-band high resolution NAIP imagery that I had put into a mosaic and am now trying to create a raster using the Mosaic to New Raster tool. When run the tool just says 'mosaic is loading'... and that it was succeeded 10 hours later. When going to the gdb that I wanted to place the raster tif into, the raster doesn't exist. It was like it never ran although ArcGIS Desktop history says it had. I had run a larger set of data through the tool on 3-band imagery and had no issue. I think I am having a similar issue to the question posed here: Mosaic to New Raster doesn't actually output anything

Does anyone know why this is happening? Is there a parameter that I am missing?


texts: ESRI Javasacript API Legend digit: Modify Layer Names

texts:
I am using the ESRI Javascript API Legend Digit in a map application and I would like to alter the display names as they appear in the legend.

I'm accessing map services with cryptic names like "Cel_clip2_V24" and I would like them to say something like "Local Farmer's Markets".

How do I control either the layer name or the legend label to achieve that?


texts:
Can anyone tell the program code or script for converting HDF files to geotiff using MRTSwath tool?

I've downloaded MODIS Files for fog estimation but i couldn't find how to process them.


texts: How to install matplotlib and pyparsing correctly on OSX-Lion?

texts:
I am trying to install the Contour plugin in QGIS 2.0.1 (I'm running OSX-Lion) and I have installed NumPy and Matplotlib from kyngchaos however when I go to install the plugin I get the following error:

  The Plugin is broken. Python said:
  matplotlib requires pyparsing

Did I misinstall matplotlib or do I simply need to find a copy of pyparsing, and if so, where is a good place to get it?

UPDATE:
I downloaded pyparsing from pyparsing.wikispaces.com and have been following the guidelines for installation. The README tell me to "Do the usual: python setup.py install" but the online instructions say to use "easy_install pyparsing". I am not sure what I am really supposed to do with this information. I have tried running the "Setup.py" file but it lacks a command but I am not sure where to add the command.

Thanks
DF