texts: Determining direction of rivers i.e. which villages lie upstream and downstream?

texts:
I want to determine which villages lie upstream and which villages lie downstream big cities built on the side of rivers in Indonesia. I have the .shp file of villages and cities, the .shp file of all rivers in the country and some .tif files showing elevation (standard SRTM rasters, although I cannot find the exact resolution right now).

I know how to identify villages located close to rivers (using spatial join with the rivers .shp file). However, I do not know how to determine the direction of the rivers. Should I discard my rivers .shp file and try to create a Digital Elevation Model from the SRTM .tif data? Using the .tif data I just did: FLOW DIRECTION + FLOW ACCUMULATION, but the computer is stuck since 40 minutes ago and even if it works I would not know how to continue.

Any suggestion?

At the moment the river polylines are a bunch of segments divided (most of the times) by nodes. In order to intersect correctly the river polylines with the village polygons I think that the river data should have one observation for each river (rather than segment): typically the river would start at the top of a mountain and would end in the sea/ocean. I think that the share of villages intersected by more than one river would be much lower than the share of villages intersected by more than one segment.

I tried to do this by dissolving the segments by one of the two node identifiers (tnode and fnode in the attribute table). Things improved, but I definitely did not manage to get only one observation for each river (I suspect the data must be also expanded somehow and not just dissolved). This is why I left the raw data in the sample and you may have had the impression that most of the villages are cut by many rivers.

I guess that I need to find a good summer school in hydrology (or an hydrology co-author) to come out of this issue.


texts:
I'm trying to clip a map of Oakland City Council Districts to the shoreline. Carto even has a nice walk through of how one could (theoretically) accomplish this, but I feel like I'm missing something basic here.

If I have a map of US States, and a map of Oakland city council districts, I want to do something like this:

SELECT oakland.*, states.the_geom
FROM oakland_city_council_districts AS oakland,
 US_states_shapes AS states
WHERE ST_Intersects(oakland.the_geom, states.the_geom)

But that particular query just gets me an error: column reference the_geom is ambiguous -- if I pull out states.the_geom I get the shapes, but they aren't cropped at all.

Realizing, with help in the answers here, that I should be using ST_Intersection (not ST_Intersects) I revised my query to:

SELECT oakland.councilmem, oakland.dist, oakland.fullname,
ST_Intersection(oakland.the_geom, states.the_geom) AS the_geom
FROM oakland_city_council_districts AS oakland,
US_states_shapes AS states WHERE oakland.dist = 3

But that query finds nothing.

So how can I use the single California polygon in my US_states_shapes table to trim the Oakland City Council Districts to the coastline?


texts: Geodetic lat/lon to ECEF cartesian: Direction to apply geoid undulation from ellipsoid surface?

texts:
I'm working on modeling the earths surface in 3D, using a combination of WGS84/EGM96 and NASA SRTM orthometric height data.

To calculate a final cartesian XYZ coordinate for a point on the earth surface, I am using the following methodology:

I convert a given (geodetic) Lat/Lon to a point on the WGS84 ellipsoid using standard math.  To convert that position to a point on the surface of the geoid, I offset along the ellipsoid normal from my ellipsoid surface point by whatever the geoid undulation is at that location according to EGM96.

In order to come to a final surface point of terrain, I then add the SRTM orthometric height to my geoid surface point, along the normal to the geoid (ie not along the ellipsoid normal!).

I'm not sure if I am offsetting from ellipsoid surface to geoid surface correctly by adding the height along the ellipsoid normal.  If you look at the following diagram I found, they seem to be adding the geoid undulation along some arbitrary direction (ie not normal to ellipsoid or geoid, just straight up).  Is the diagram wrong, or am I wrong to add the undulation height along ellipsoid normal?

Here is the image I am referring to, and as you can see, the refer to the geoid undulation offset as "N" in the diagram:

Source: http://kartoweb.itc.nl/geometrics/Reference%20surfaces/refsurf.html

Thanks for any help in advance!


texts: Getting distances from points to polygons using QGIS?

texts:
I'm looking at finding the distance in QGIS between a few things.

Basically, I'll have a points layer, and a polygon layer. Both will have several features in them.

I need to find the distance from each point to each polygon. I can potentially get the centroid and use that (and I think that's what will probably be used), but bonus points for the nearest edge of each polygon to the point in question.

Also, for some reason, v.distance doesn't seem to be appearing in my processing toolbox.

Also tried getting centroids and using the distance matrix, however that's giving me incorrect values.

EDIT - I've done a quick trial (still waiting for the actual data) using the standard distance matrix and measuring to centroids. This works well, but unfortunately the measurements returned are slightly out. This will probably be down to the coordinate system, and I'll have to see what the actual data is like when I get it.

Unfortunately, MMQGIS's hub distance only finds the nearest (and even then, not to the edge of the polygon but to a point within, which doesn't appear to be the centroid).

In any case, the people after this data say that the centroid is okay - they just need to get me the data (which is a tad frustrating, as it is due today).


texts: PostGIS: Ensure Intersection of 2 Polygons results in a Polygon

texts:
I want to ensure that the result of a ST_Intersection(polygon1, polygon2) is also a polygon. I initially started by filtering with a ST_Overlaps(polygon1, polygon2) but that seems to only ensure that both features intersect while neither feature contains the other.

Besides filtering by checking the result of the intersection operation e.g. ST_Dimension(ST_Intersection(polygon1, polygon2)) > 1 is there a way to filter polygons by ensuring that they share some internal space that is not a boundary?

ST_Crosses(polygon1, polygon2) actually does perform the filtering that I want.


texts: Aligning dataframes to exactly match information in another dataframe when sizes of 2 dataframes are different?

texts: I have a dataframe which is using data-driven attribute (say zone). Now this layer has grids with specific numbering but since it is data driven, each time the zone changes, the grid numbering will change as well. To make it easier for the user to see the grid numbers, i created another dataframe with just the grid numbers and tried to put it behind my original map. But, i am not able to make them align with one another perfectly such that the grids match between both the dataframes. I have tried setting the extent of the background dataframe to be based on the original dataframe's extent (100%) but it doesn't help as the sizes of 2 dataframes don't match. I am creating a mapbook with 100's of maps and would like to set it to match the grids automatically for each zone.  I am attaching the sample image for easy reference of the question.


texts: Automatically convert from latlon to planar

texts:
I work with library(sp) and points in latitude longitude, with this projection:

proj4string(pts) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

But some functions (such as sp::spample) require planar coordinates (for good reasons).

Sometimes I work with points in Europe, sometimes in Asia, etc.

Is there a best practice to project to planar coordinates without having to manually specify in which epsg-zone the coordinates are?

This will do it (but knowledge that points are in Finland must decide "2393").

pts.transfor <- spTransform(pts, CRS("+init=epsg:2393"))

How can "2393" be found automatically in R?


texts: Aggregating rates represented by points within X radius using QGIS

texts:
the point represents rates like 10%, etc. (this was also pretty clear from 2 years ago) -- been trying to do this for 2 yrs...

the points dont represent numbers of events since they are rates (Finding largest number of observations within a given radius)

but the points have geo ids so they are attached to location, not distance (whatever distance is to imply) -- update: it represents z-value but we're just dealing with rates


texts: Stop layer from showing up in legend on composer

texts:
When I add a legend in a map composer in QGIS, all layers that are defined in the project will be shown in the legend, no matter if they are set to be displayed or not, I also have some background maps I would like to never show up in the legend.

So:
1) Is there a way to tell Qgis not to list layers that are not displayed in the legend?
2) Is there any way to tell Qgis to never list a layer in a legend?

(In my case, I would like to never ever show the background layers or the non-displayed layers the legend)

(using Qgis 3, but have had the same issue for quite a few 2.x versions as well)


texts: How to calculate pixel area of a clipped image in Google Earth Engine?

texts:
I calculated NDVI using Sentinel 2 image collection by mapping a function over the image collection and then returning the expression inside the function as a band .addBands(expression).rename('NDVI')
Upon printing, every image had a new band 'NDVI' in the console.
Next I created a mean of image collection to make a single image and selected the 'NDVI' band I just created.
Next I clipped the NDVI image using a shape file.
Now all I want is to print the pixel area of the clipped NDVI image inside Google Earth Engine console.


texts: I am working on remorte Desktop where arcGIS 9.3.1 is installed. It was working perfectly for Long time. Currently it shows a funny error like: you do not have enough memmory to run the application and stop before opening arcmap or sometimes stop in the middle of a Tool progress...


texts:
I want to use my custom created mapbox basemap in OpenLayers application,  I have added it at TileLayer with XYZ source, but the quality is getting blurred

This is the image from application

This is how it looks in Mapbox Studio.
My code
new TileLayer({
          source: new XYZ({
            url:  'https://api.mapbox.com/styles/v1/johnculturemee/cknzq8w4c4yq617jgv5kwg9tm/tiles/{z}/{x}/{y}?access_token=<my_key>',
          }),
        })

What might be the reason?


texts: Selecting features from multiple layers for export QGIS 3.6

texts:
I have multiple features I want to select over multiple layers and then export these selected features to a new layer.

The 'Multiple Layer Selection' plugin is not supported by the QGIS 3.6 so that is not an option.

I know I can create one massive layer file and use expressions to display and then select only the features I want but this is a cumbersome process compared to being able to use the default (rectangle) select tool across multiply layers.

Is there potentially a way to use the identify feature tool to create an attribute table that contains all the identified features across all layers? That would add an extra step to the selection process but is still better than using multiple expressions on a massive layer file.


texts: Setting Python Toolbox Parameter to None without making current selection value unselectable thereafter?

texts:
I'm making a python toolbox and it has an input parameter of type GPFeatureLayer. Let's call it my_param.

In updateParameters when my_param.altered, its value is added to a value table. Then I'd like to set my_param.value = None, so that the selection is blanked out and same value could be immediately selected again.

As it stands right now, setting my_param.value = None means that that layer cannot be chosen again.

Why is this happening, and how can I achieve what I'm trying to do?

Here is some example code for you to use:

def getParameterInfo(self):
    my_param = Parameter(displayName='My Param',
                         name='my_param',
                         datatype='GPFeatureLayer',
                         parameterType='Optional',
                         direction='Input')
    return [my_param]

def updateParameters(self, parameters):
    my_param = parameters[0]
    if my_param.altered:
        (function that appends my_param.value to another table, not important for this problem)
        my_param.value = None
    return


texts: GDAL gdalwarp usage on QGis - how to set cell position of the clipped raster to be equal to the origin of the source cell

texts:
I work on a Windows QGis 2.14 platform. I need to clip a raster with a vectorial polygon. When I use gdalwarp I cannot find where is the setting for establishing the output graticule position respect the source raster. I know that in the ESRI ArcGIS clip raster tool is there a command like "maintain the same graticule as the source". Without this option gdalwarp gives me back a raster with cell shifted (offset) respect the source. In fact the new origin of the output raster is determinate by the rectangular envelope of the polygon (obviously).

Suggestions for me, also workaround if necessary?

Similar questions didn't answer this question.


texts: How to create group of coordinates and check if a coordinate is within the group?

texts:
I am looking for a way to relate two or more types of open data from the Valencia Open data portal. This portal offer data like a traffic jam, noise, and air pollution. All these data have coordinates based on the locality of real sensors. So, I am using GeoTools which is a java library, to work with these data.
Now I want to create some polygon groups inside Valencia city in order to say if some coordinates are close to each other. For example, I imagine that if a coordinate XY and XY' are in the same group, it is likely that they are close to each other.

How can I create these groups? I am searching here but I don't have too many ideas on how to proceed on this.

For example, let's say I have 5 points and I want to create sets with points that are close to each other in at least 3 meters. My result would be 2 sets like below:

point A: 44, 33
point B: 44, 34
point C: 45, 34
point D: 54, 53
point E: 55, 55

set 1: points {A, B, C}
set 2: points {D, E}


texts: Creating vice versa multipoint feature within ArcGIS for Desktop?

texts:
The help says about a multipoint feature:

  Sometimes, you need to create a feature that has more than one physical part but only references one set of attributes in the database. These are called multipoint features for point or multipart features for line and polygon features.

I have an reversed situation. I need to create a single feature with multilined attribute table.

For example I have a building under construction. Every day it is raised at specific height. The contour of this building is always the same, but the height is changing every time.

So I need few attribute lines to describe such building.


texts: ArcGIS: How can frequently used colors be included in the 8-bit palette when exporting to TIFF to avoid dithering of large parts of the image?

texts:
ArcGIS Desktop 10.7.1 can export a layout to TIFF with 8-bit palette. But the 256 slots in the palette are filled with gray-scale colors only used for one or a few pixels when black text on white background is anti-aliased. Then there are no free slots in the palette for the important colors used in large parts of the image. Therefore, those are dithered. I can see 4 solutions but do not know how to achieve any of them:

Tell ArcGIS to be smart when creating the palette by filling the slots with colors in order of the number of pixels that should have that color. The color that is most common will be added first, then the second most common and so on.
Tell ArcGIS to not anti-alias text when exporting to a raster graphics format.
Manually tell ArcGIS that certain colors should be in the palette and let ArcGIS fill the rest of the slots.
Make ArcGIS use a manually crafted complete custom palette where the important colors are included, and the remaining slots are filled with colors selected to be suitable for anti-aliasing of the text in the layout.


texts: Leaflet VectorGrid interactivity doesn't work if Leaflet.markercluster is used

texts:
This is a cross-post from SO. I am trying to use Leaflet VectorGrid interactivity option for click and mouse events:
const vectorGrid = L.vectorGrid.slicer(geoJsonDocument, {
    rendererFactory: L.canvas.tile,
    vectorTileLayerStyles: {
      sliced: geoJSONStyle(false)
    },
    maxZoom: 22,
    indexMaxZoom: 5, // max zoom in the initial tile index
    interactive: true
  });

  vectorGrid.on("mouseover", function (e) {
    console.log("mouseover");
  });

  vectorGrid.on("click", function (e) {
    console.log("click");
  });

However, interactivity doesn't work if Leaflet.markercluster is used.
I created a codesandbox.
If you comment map.addLayer(mcg);, interactivity works.


texts:
I'm trying to use both PostGIS and GeoTools with the same ESPG:3857 projection, however, it looks like the projections in PostGIS and GeoTools differ.

This is what I have in PostGIS:

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]

Here's GeoTools:

PROJCS["WGS 84 / Pseudo-Mercator",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic latitude", NORTH],
    AXIS["Geodetic longitude", EAST],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Popular Visualisation Pseudo Mercator", AUTHORITY["EPSG","1024"]],
  PARAMETER["semi_minor", 6378137.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["central_meridian", 0.0],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH],
  AUTHORITY["EPSG","3857"]]

I've tried using the PostGIS one with GeoTools, but it cannot parse it due to the EXTENSION element. Of course, I could update PostGIS with the GeoTools one, but can anyone tell why they are different?

I'm using GeoTools 10.0, Postgres 9.3 with PostGIS 2.1. I'm just trying to convert a EPSG:4326 point into EPSG:3857, and vice versa. The code works:

    CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
    CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:3857");
    MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, false);
    GeometryFactory geometryFactory = new GeometryFactory(new PrecisionModel(), 4326);
    Point point = geometryFactory.createPoint(new Coordinate(lon, lat));
    Point targetPoint = (Point) JTS.transform(point, transform);

but the targetPoint that I get is different from ST_Transform(point, 3857) that I get from PostGIS.

So for example, for

(41.937832, -87.642970)

I get from GeoTools with code above:

POINT (4668498.103213854 -24771604.62292875)

and from PostGIS

select ST_AsText(ST_Transform(ST_SetSRID(ST_Point(-87.642970, 41.937832), 4326), 3857))

POINT(-9756370.79201015 5151671.52336743)

Any solution on how to make these two work together? I'd really like to convert between projections in both SQL and Java.

Thanks!


texts: rendering style symbol (star) separately, not on the map, Openlayers

texts:
I try to create a legend for the map separately in a different div.

How could I render a symbol (star e.g.) somewhere else in html, not only on the map with OpenLayers?

I guess I have to call OpenLayers.Renderer.symbol.star and display it elsewhere, but I couldn't figure out which method should be used...

The most perfect solution would be to get the exact styling properties like on the map...
Is there a way for it?


texts: Apply raster data to points from multiple years

texts:
I am quite new to QGIS so please bear with me.
My data set is of presences data with latitude, longitude, month, and year of cetaceans. The time frame is from different months and from the year 1990 to 2018. I have multiple raster files for sea surface temperature for each month and year and I want to apply the sea surface temperature with the appropriate data point.
I have applied bathymetry data and slope to my points by using the following code in the field calculator
raster_value('rast', 1, make_point($x, $y))

Is there any way I could do change or add to this code to apply the data to the correct year and month?
or
Does any one know of a simpler way of doing this?


texts:
When I write

x <- c('2007-05-05', '2006-12-31')
strptime(x, "%Y-%m-%d")$year

I get the result

[1] 107 106

I thought that POSIXlt started at 1970, but this appears to start at 1900?


texts: Changing scale bar interval and other settings on individual Data Driven Pages?

texts:
Sometimes I'll have a MXD set up as a mapbook using Data Driven Pages, but I'd like to display certain pages at a different scale than others in the same mapbook. I can do that using the Data Driven Scale option in Data Driven Pages, but when those pages are viewed or rendered, the scale bar uses the same units and divisions (for example, two 1000-foot sections) for all pages and is simply stretched or compressed to correctly match the display scale of each page. So if I have a mapbook where most pages are 1:12000, and one page is 1:6000, that one page will have a scale bar that is twice as long as the other pages, because ArcGIS is using the same settings (division units, division value) for the scale bar on every page.

Is there a way to change this, so that I can set the scale bar parameters independently for each page?

For example, if I go from a scale of 1:24000 to 1:10000, I might want to use increments of 0.25 miles and 500 feet respectively.

Currently what I'll do if I have, let's say, an 80-page mapbook with a few pages at a different scale than the rest, I'll render the majority of the pages first, then change the scale bar settings manually for the few pages that would look better with a different scale bar. But if someone comes along and opens the MXD after me and tries to render the same mapbook, they may not take that step and the scale bars will not be right.


texts:
We are currently migrating to ArcGIS Pro 2.2, with new add-in to improve our workflow.
Our Add-in have been shared on a common network drive so far. But our main concern is that we are not admin on our machines, so when we need to install an add-in the IT guy needs to come.

A way to share any add-in and to update them would be a silent install. I found this thread on Geonet for ArcGIS desktop.
The way described is to change the default path to the add-in to a new path with the proper registry key:

  \\HKEY_LOCAL_MACHINE\SOFTWARE\Wow6432Node\ESRI\Desktop10.2\Settings\AddInFolders and set the value

Is it possible to use the same method for ArcGIS Pro add-in?


texts: Why are there two Draw Effects buttons in Layer Styling Panel in QGIS?

texts:
At 3.10.3, I notice that the Draw Effects button is located twice in the Layer Styling window, as shown in this screenshot:

I've spent time testing them, and they both work equally well; that is, they seem to do the exact same things.  One button's functionality is no different than the other.  Additionally, if both are selected, they work together additively (however, the additive effects of the second button could be replicated with a single button).
So I'm curious why there are two buttons. The docs don't discuss this duplication.
In summary, what is the purpose of a second Draw Effects button?
EDIT: I also have access to 2.18 which has only one Draw Effects button, within the Layer Rendering section.


texts:
I have a control that draw line. What I want is that once the line is drawn the control should be deactivated.

new GeoExt.Action(Ext.apply({
    text: 'Draw Line',
    control: new OpenLayers.Control.DrawFeature(
        this.layer,OpenLayers.Handler.Path, {
        eventListeners: {
            featureadded: this.myFunction,
            scope: this
        }
    })
}, actionDefaults))

myFunction: function(evt) {//Deactivate Control Here }


texts: How does ArcGIS Desktop recognize WMS CRS?

texts:
I'm using ArcGIS Desktop 10 and working on a WMS service with Mars data. Currently the GetCapabilities shows <SRS>EPSG:4326</SRS>. I assume ArcGIS 10 reads this and therefore the WMS has GCS_WGS_1984 as a Geographic Coordinate System?

I would like to use a custom URL definition because I am working with data from Mars. Is it for example possible to add a URL with GML such as http://spatialreference.org/ref/iau2000/49900/gml/ ?


texts: QGIS set fixed ID name with autoincrement in data attribute table

texts:
I am going to create a bulk of lines in QGIS.
One line will have allocated ID, which looks as follows:
DT.C.10_Strand
Thereafter each "Strand" will have the number of order like you see below:

After adding a new row I need to increment my ID number by one.
How can I define it in advance? Is there any solution in i.e. field calculator or somewhere?
So far I have to write down every time the ID name. Can I have it ready and only change the number? or increment id as we can do in Excel?


texts: Adding geoprocessing results as layer to dataframe using ArcPy?

texts:
I am still new to ArcPy.

I'm running a really simple script, it selects a feature from a dataset (one town out of many), saves the data, and adds it to the (empty) dataframe.

    map = mapping.MapDocument("CURRENT")
    df = map.activeDataFrame
    where_clause = "TOWN = " + "'"+GetParameterAsText(0)+"'"

    arcpy.MakeFeatureLayer_management(statewide_towns, "TOWNS")
    arcpy.SelectLayerByAttribute_management("TOWNS","NEW_SELECTION", where_clause)
    arcpy.CopyFeatures_management("TOWNS", "OUTPUT_Single_TOWN")
    arcpy.MakeFeatureLayer_management("OUTPUT_Single_TOWN", "TOWN_BND")
    mapping.AddLayer(df,"TOWN_BND","AUTO_ARRANGE")

Throws me this error:

Traceback (most recent call last):
  File "C:\path\to\myTool.py", line 58, in <module>
    mapping.AddLayer(df,"TOWN_BND","AUTO_ARRANGE")
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\mapping.py", line 49, in AddLayer
    assert isinstance(add_layer, Layer)
AssertionError

Failed to execute (myTool).


texts: ArcGIS: Get coverage of Polygon A in Polygon B

texts:
We have dissemination area (DA) units for a city [Polygon B], and shapefiles of all parks in the city [Polygon A].Polygon A overlaps Polygon B in some areas.

We initially wanted to calculate the percentage park cover in each of the dissemination areas but now we'd just be happy to create a feature class for each of the dissemination areas which includes a field in the attribute table indicating the total cover of park within its boundary.

We've tried pretty much all of the overlay tools listed in this webpage:
http://desktop.arcgis.com/en/arcmap/10.3/analyze/commonly-used-tools/overlay-analysis.html ,but they haven't been able to give us a feature class with the same number of DA's as Polygon B, with the park data nested inside.

Is there a glaringly obvious tool we've missed?

It seems to be a pretty easy task but we've really hit a bit of a wall with it.

Here is a screenshot. Polygon A in green, Polygon B represented by the red lines.


texts: Error using 'Join attributes by location' via Python

texts:
Trying to run a simple script using qgis:joinbylocation but getting either errors or it's not processing correctly. I don't want to take a summary of the attributes, just the attributes of the first located feature and I want to keep all records. The parameter selections should be correct but it's the sum,mean,min,max,median that's causing problems. I tried removing them but python wasn't too happy. Also (for some reason) python doesn't recognise mean or median but is fine with the other statistical terms.

Any advice please?

##Test=name
##File_1=file
##File_2=file
##Output=output vector

mean=None
median=None

processing.runalg("qgis:joinbylocation", File_1, File_2, 0, sum,mean,min,max,median, 0, 1, Output)


texts: Exporting from PostGis using Ogr2Ogr - "Unable to open datasource"

texts:
Wrestling with Ogr2Ogr again here. I have managed to use it to import a geojson file into a PostGis database, can query and manipulate the database layer. However, when I try and export the database table to a new file I get the message: "Unable to open datasource 'myexport.geojson' with the following drivers"

GeoJson is in the list of drivers, and the layer is in PostGis.

  .\ogr2ogr -f "GeoJSON" "PG:host=**** port=**** dbname=**** user=**** password=****" "C:\Temp\output.geojson"


texts:
I've inherited some scripts that (to my mind) "abuse" geocoding to position points based on a LocationID field.  As the name implies, this field is a unique ID that really doesn't need to be "fuzzy matched"; either it matches or it doesn't.  Besides that, using the geocoding tools for this process is painfully slow; some 5 minutes or more to geocode less than 100k records.

I've tried using joins, with not much luck:

Joining the tables to the reference points with AddJoin: with "KEEP_ALL" results in all the points, some with no data; with "KEEP_COMMON" results in only matching points, but loses data where multiple rows have the same LocationID.
Joining the points to the tables with AddJoin results in tables, not point features.
Joining either way with JoinField is even slower than geocoding, and has the same issues as above.
None of these methods produces a new feature class as geocoding does; that takes a CopyFeatures either before or after the join.  (SpatialJoin does, but it's out of the question since one input is a table....)

So, is there any faster way to do this using GP tools – if not, that might explain why it is the way it is – or should I look to a lower-level solution (arcpy.da or the like)?

(To clarify, the output I'm looking for is a point feature class:  each point has data from a row in the input table, but there may be multiple points in a given location.)


texts: Using flat (no data) areas of aspect-raster in QGIS raster calculator?

texts:
I´m playing around with the Raster tools in QGIS, following more or less the qgis-tutorial but using my own data. When I calculate the aspect of a DEM, flat areas are (logically) left without data. However, when I try to make an addition with another rather afterwards in the raster calculator, the flat areas (no data in my aspect layer) are never taken into consideration, even if the raster I use for addition has data in this area.

For example: Let´s say I´m looking for places facing south (135-225°) with slopes less than 15%. To find these places I have two rasters:

Raster A: aspect 135-225° has the value 1, the rest 0 and flat areas
without aspect have no data (the appear transparent)
Raster B: slopes less 15% have the value 1, the rest 0

If I combine the two rasters in the Raster Calculator (e.g. like this: "ANP_aspect_sur@1" = 1 AND "ANP_slope_15@1" = 1), the output file has no data where Raster A has no data, independent of Raster B.

What am I doing wrong?

I thought about setting flat areas to 0 (-zero_for_flat). But in this case 0 would stand both for north in regards of the aspect and for nodata.

As far as I know ArcMap assigns -1 to flat areas.


texts:
I am trying to draw a polygon dynamically.I get the latitude and longitude from two different entities.Hence I concatenate them as below

var greenArr=new Array();
var point = new esri.geometry.Point(-85.91,25.68,this.parent.esriMap.spatialReference);
var greenPolygon = new esri.geometry.Polygon(this.parent.esriMap.spatialReference);
if (type == 'Green')
{
var coords = '[' + parseFloat(longitude)+ ','+ parseFloat(latitude)+ ']';
greenArr += coords + ",";
}
greenArr="["+ greenArr +"]";
greenPolygon.addRing(greenArr);
var result=greenPolygon.contains(point);

Here the contains method is not working.it returns false always


texts:
As I'm using multiple functions in one process, I would like to pass results of each function to the next one using an In-memory File. I'm using Rasterio and I saw that there is a way to create In-memory files in the documentation, but unfortunately it seems like it's not working like usual files. I cannot use open(temporaryfile).meta.copy() to get infos (GTIFF in my case) of the file.

Is there a way to correctly store everything as in a Gtiff file but in a memory file or do I have to process in another way ?


texts: Add raster paths to shapefile with Python (ArcGIS)?

texts:
I have a problem and maybe someone is able to help me or give me some hints how to start.

Situation:
I have a lot of raster and of every raster I have the center coordinate. Out of this coordinate I build a shapefile to have the information about the approximate location of each raster. Now I included the path for each raster in the attribute table of the shapefile with the field calculator.

Problem:
Some of the raster are exactly georeferenced with a word-file. I want to add a new attribute field to show if there is a georeference or not. For this I can use the following code (at the moment a input cursor is missing):

from arcpy import env

arcpy.env.overwriteOutput = 1

env.workspace = r'Z:\path'
path = env.workspace

# Create a table and add fields
table_name = "TempTable.dbf"
arcpy.CreateTable_management(path, table_name)
arcpy.AddField_management(path + "\\" + table_name, "RasterName", "TEXT")
arcpy.AddField_management(path + "\\" + table_name, "cs", "TEXT")

# Get the list of raster files
rasters = arcpy.ListRasters()

# Loop through the files and perform the processing
for ras in rasters:

    # Determine if the input has a defined coordinate system
    dsc = arcpy.GetRasterProperties_management(ras, "TOP")

    if float(str(dsc.getOutput(0)).replace(",", ".")) > 100:

        row.cs = "y"     # Adds "x" to the row where the image has a georeference

    else:

        row.cs = "n"     # Adds "n" to the row where the image has no georeference

Now I do not know how to add the information in the right way to the attribute table (maybe a join?). I hope there is someone who is able to help me with this.


texts: Why does adding DWG from AutoCAD Civil 3D 2012 freeze ArcMap 10.1?

texts:
I use AutoCAD Civil 3D 2012 and ArcMap 10.1.

I am trying to add a .dwg file into ArcMap and no matter how I try (drag and drop from Catalog window, Add Data) ArcMap freezes and locks up.

Any ideas?


texts: Visualization of a big graph in 3D from postgresql/postgis

texts:
We have a graph in 3D in a PostgreSQL/PostGIS database. It contains something like 50'000 edges and 50'000 vertices, with x-y coordinates and floors.

We would like to visualize this graph with straight lines as edges in 3D and to add on top of this few more time-dependent data related to the graph.

What would be the best technology to do this? Anyone with similar experience?
We would prefer this technology to be:

open source
as easy as possible
but still able to manage this graph and time-dependent data

Would Processing be able to manage this?
Or would it be better to go for OpenGL?
Is there any other technology that would be more appropriate for this?
And how long would it take, if you've any experience, to visualize the graph with each technology?

I've hear that processing already reduces the frame rate with 10'000 links. Also, it seems interesting to use a bit of both. OpenGL seems nice since you can program it in any language (c, c++, objective-c, c#, java, etc.) and results could already be obtained in one week.


texts: GeoWebCache spending so much more time than WMS service to create tiles

texts:
I'm trying to undestand what is happening with GeoServer/GeoWebCache integration.

Anyone has experienced GeoWebCache spending so much more time than WMS request to generate tile ?

In a simple test in my environment using embedded GWC in GeoServer, I also tested GWC in standalone mode, when I call a tile in GeoServer direct using WMS service the tile is returned in 18s, this same tile passing by GWC is spending more than 1min and return an error "400: ImageIO.read(InputStream) returned null. Unable to read image" because my WMS limit is 1min. In this specific case the tile has average of30.000 features. In other cases with less features GWC works well and return tiles without problems.

I thought that GWC just call GeoServer WMS service save the image returned and return the same to the client, so in this case if GeoServer spent 17s to return image to GWC, why the GWC can´t process this image in less than a minute ?

Something strange is that calling WMS Service directly the image is returned in 17s but using GWC the same request exceeds the WMS time limit that is 1min in this case and a error is throw in GeoServer log.

I'm using GeoServer 2.9.2 with GWC 1.9.1 and has tested GWC 1.9.2 standalone with the same results.

Anyone have any idea about it ?


texts: Accuracy of measurements in Google Earth

texts: Can anyone point to a resource that lists the accuracy of measurements taken in Google Earth?


texts:
I have a web mapping application built using the Web App Builder function which contains a print widget. I am using the default ESRI print service (ie the URL is blank in My Organisation - Edit Settings - Utility Services - Printing - URL is blank). The feature layers (created in my ArcGIS Online account) are all publicly shared. I get en arror when I go to use the print function to create an image or pdf.
This is the URL to the web mapping application http://missgis.maps.arcgis.com/apps/webappviewer/index.html?id=1121af9fcf4f47399ebb0873993689ff

Any ideas?


texts: Display basic geoJSON polygon on OpenLayers map

texts:
For some reason my basic polygon, which I've verified on geojsonlint.com, won't display correctly on the map. Looks like a projection issue, but I can't determine what I missed. I'll paste the code below, and here is the jsFiddle:

http://jsfiddle.net/kzex47za/

    var geoJSON = {
  'type': 'Feature',
  'properties': {
    'id': 'shape'
  },
  'geometry': {
    'type': 'Polygon',
    'coordinates': [[
      [-95.0144,39.0061],
      [-95.0144,38.9842],
      [-94.9862,38.9842],
      [-94.9862,39.0061],
      [-95.0144,39.0061]
    ]]
  }
 }

var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM()
        }),
        new ol.layer.Vector({
          source: new ol.source.Vector({
            features: (new ol.format.GeoJSON({
              featureProjection: 'EPSG:3857'
            })).readFeatures(geoJSON)
          }),
          style: new ol.style.Style({
            stroke: new ol.style.Stroke({
              color: 'rgba(29, 130, 29, 1)',
              width: 2
            }),
            fill: new ol.style.Fill({
              color: 'rgba(0, 255, 0, 0.45)'
            })
          })
        })
    ],
    view: new ol.View({
        center: ol.proj.fromLonLat([0,0]),
        zoom: 1,
        projection: 'EPSG:3857'
    })
});


texts:
I am trying to change the CRS of a web mapping service in Norway. The projection is EPSG:25832.

Here is the WMTS service: https://kartkatalog.geonorge.no/metadata/norge-i-bilder/norge-i-bilder-wmts-euref89-utm32/61ed45e1-2ef3-4f94-873c-e9dcb6be4aa6

Here is the example from Proj4Leaflet Documentation:

var crs = new L.Proj.CRS('EPSG:2400',
  '+lon_0=15.808277777799999 +lat_0=0.0 +k=1.0 +x_0=1500000.0 ' +
  '+y_0=0.0 +proj=tmerc +ellps=bessel +units=m ' +
  '+towgs84=414.1,41.3,603.1,-0.855,2.141,-7.023,0 +no_defs',
 {
  resolutions: [8192, 4096, 2048] // 3 example zoom level resolutions
 }
);

The WMS loads using the code below when the CRS information is not being used. The map centers on northern Norway but it should be centering near Oslo.

There is something wrong with the CRS information then or how I am using proj4.js?

var crs = new L.Proj.CRS('EPSG:25832', '+proj=utm +zone=32 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs',
{
  resolutions: [
     131073, 65537, 32769, 16385, 8193, 4097, 2049, 1025, 513, 257,
129,
65, 33, 17, 9, 5, 3, 2
  ],
  origin: [-2000000.0, 9045984.0],
  bounds:  L.bounds( [-2000000.0, 3500000.0], [3545984.0, 9045984.0])
})

var key ='****************'

var map = L.map('mapid',{
zoomControl: false,
zoomSnap: 0.1,
// crs: crs,
continuousWorld: true,
worldCopyJump: false
}).setView([59.877812, 8.590628], 5);

var gkt = '***************';
var url =
'http://gatekeeper{s}.geonorge.no/BaatGatekeeper/gk/gk.nib_utm32_wmts_v2?&gkt='
+ gkt
+'&layer=Nibcache_UTM32_EUREF89&style=default&tilematrixset=default028mm&Service=WMTS&Request=GetTile&Version=1.0.0&Format=image%2Fpng&TileMatrix={z}&TileCol={x}&TileRow={y}';

var wms = new L.TileLayer(url, {
   subdomains: ['1', '2', '3'],
   maxZoom: 16,
   minZoom: 0,
   attribution: ''
 }).addTo(map)


texts: How to build a website or mobile app with a routing capability?

texts: I want to build a website/mobile app for my school that shows the direction between two points  similar to Google's directions API, but not with Google maps. I've got map shapefiles and GPS  coordinates of buildings in the school and some programming  knowledge (javascript, php,  mysql), Please how do i go about it, I'm quite new to GIS  programming.


texts: Calculate a river system with the field calculator in QGIS

texts:
I would like to create a river system with water flow based on the following attribute table:

The river system is divided in different segments (the ID to every Segment is called "SegmentNo").
The column "Downstream" shows the ID from the next downstream segment.
Not every segment have a water flow value. Only these segments with a waste water treatment plant (WWTP) in his environment.
Now I would like to compute the water flow for each segment in downstream direction with the field calculator.
Unfortunately I am not a hero in Python, so I can't create new functions.
Probably it isn't necessary too.
I have found a "When Then" solution which add up the water flow. But not in combination with the downstream segments. Because not every segment have only one neighbour.
Can I solve the problem with the given functions?


texts: script or code to create map package arc gis , mpk

texts:
I am looking for some help or if someone has a script or tool to loop through a given folder of MXD's and Create MPK files with the same name.

i don t have so much experience of creating scripts , or using the python window. I found the previous question below

enter link description here

Where it details the ESRI , help page for data managment

enter link description here

I tried inputting the python window code changing the code for my examples but i keep getting different error message.
I am not sure what i am doign wrong. I m not sure about setting up the parameters for an actual tool . i was wondering if maybe someone had already done this.


texts: Reformatting Dictionary Object to rows and columns using Google Earth Engine?

texts:
I'm new to javascript & GEE, trying to reduce a global FeatureCollection by computing sums of two specific properties, grouped by country (using this helpful link). The output is a Dictionary Object that contains all the info I need.

I'd like to reformat the data into rows and columns and then export it as a CSV.

What is the simplest way to do this?

// Call up WDPA data, filter

var dataset = ee.FeatureCollection('WCMC/WDPA/current/polygons')
              .filter(ee.Filter.or(
              ee.Filter.and(ee.Filter.eq('MARINE','0'),
                            ee.Filter.eq('IUCN_CAT','Ia')),
              ee.Filter.and(ee.Filter.eq('MARINE','0'),
                            ee.Filter.eq('IUCN_CAT','Ib'))));

// Compute sums of protected area, grouped by country

var sums = dataset
  .filter(ee.Filter.and(
    ee.Filter.neq('GIS_AREA', null),
    ee.Filter.neq('REP_AREA', null)))
  .reduceColumns({
    selectors: ['REP_AREA', 'GIS_AREA', 'PARENT_ISO'],
    reducer: ee.Reducer.sum().repeat(2).group({
      groupField: 2,
      groupName: 'country',
    })
});

print(sums);

I can export as a CSV from there, but it obviously comes out in Dictionary form (i.e. not terribly useful in a spreadsheet). This post almost addresses the question, but their initial output is a FeatureCollection, and mine is an Object. I couldn't get the chosen answer to work for me, and I'm just going in circles.


texts: Time series of NDVI calculation from Sentinal in earth engine

texts:
I am new in google earth engine, i am trying to plot the time series NDVI value and images as well, i found a script and i tried to modified but its not working properly. Below is the code

its giving error that, Error generating chart: The image collection is empty.

I just changed here the var geometry, and i imported the KML file using fusion table and paste the link here. How to display the statistics graph and images also for selected region

//Import GEE Feature Collection (My AOI kml)
    var geometry = ee.FeatureCollection('ft:19HmlYA6D2W-soGuCJGOvdEutyTNdZFK1tUFz2JAM');

// Create image collection of S-2 imagery for the perdiod 2016-2019
var S2 = ee.ImageCollection('COPERNICUS/S2')

//filter start and end date
.filterDate('2016-02-01', '2019-02-28')

//filter according to drawn boundary
.filterBounds(geometry);

// Function to mask cloud from built-in quality band
// information on cloud
var maskcloud1 = function(image) {
var QA60 = image.select(['QA60']);
return image.updateMask(QA60.lt(1));
};

// Function to calculate and add an NDVI band
var addNDVI = function(image) {
return image.addBands(image.normalizedDifference(['B8', 'B4']));
};

// Add NDVI band to image collection
var S2 = S2.map(addNDVI);
// Extract NDVI band and create NDVI median composite image
var NDVI = S2.select(['nd']);
var NDVImed = NDVI.median(); //I just changed the name of this variable ;)

// Create palettes for display of NDVI
var ndvi_pal = ['#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b',
'#a6d96a'];

// Create a time series chart.
var plotNDVI = ui.Chart.image.seriesByRegion(S2, geometry,ee.Reducer.mean(),
'nd',500,'system:time_start', 'system:index')
              .setChartType('LineChart').setOptions({
                title: 'NDVI short-term time series',
                hAxis: {title: 'Date'},
                vAxis: {title: 'NDVI'}
});

// Display.
print(plotNDVI);

// Display NDVI results on map
Map.addLayer(NDVImed.clip(geometry), {min:-0.5, max:0.9, palette: ndvi_pal}, 'NDVI');


texts: Treatment of Landsat 7 SLC error in Google Earth Engine?

texts:
Google Earth Engine distributes Landsat 7 image composites of indeces such as NDWI and NDVI that, when averaged on a shapefile, give back a coherent value.

But how is SLC failure treated by it https://landsat.usgs.gov/slc-products-background?


texts: Uploading MDB into postgis with non geometry tables in the mdb

texts:
In osgeo4w I am uploading an mdb into postgis. It works perfectly with all the geometry content within the personal geodatabase, but it does not take the tables from the personal geodatabase.

Is there any way of adding some functionality to the following osgeo4w command in order to also extract the non geometry tables from the personal geodatabase?

ogr2ogr -overwrite -f "PostgreSQL" PG:"host=... dbname=... user=... password=..." filename.mdb

UPDATE

When using ogrinfo it does not recognize tables without a geometry. They are in the .mdb when you open it in Arcmap.

ogrinfo filename.mdb


texts: Finding and changing measurement units using QGIS?

texts:
Using QGIS 2.4, I want to create 1,000m buffer around an object, but when I enter 1,000 as the distance I get something that looks more like a 100,000m buffer.

I can enter 0.01 to get what looks like a 1,000m buffer, but I would rather work in single metres.

How do I change the settings so that the number I enter is in single metres?


texts: Creating hole in polygon while using 'auto complete polygon' in ArcMap?

texts:
In ArcMap, it is possible to create a hole in a polygon by first drawing the outer limits of the polygon, clicking 'finish part' and then creating a hole somewhere in the middle of the polygon.

It is also possible to create adjacent polygons with the Auto-Complete Polygon by clicking somewhere within an existing polygon, drawing the outer limits of the new polygon and then returning to the existing polygon.

But is it possible to combine these methods, so that after returning to the existing polygon with the Auto-Complete Polygon tool you hit 'finish part' and then move on to create a hole in the middle of the new polygon?

When I try this, after having completed all of the steps above (I have drawn the outer limits of the new polygon, returned to the existing polygon, clicked 'finish part' and finally drawn a hole inside of the new polygon) I end up with only the outer limits of the new polygon. After double-clicking or hitting f2 to complete the polygon the hole has disappeared.


texts:
I am trying to make a simulation using HEC-GeoHMS in ArcGIS and I have encountered a problem near the final steps. In the drop-down menu HMS, when I click "Check Data", a report is being done and displayed in Notepad with the layers and information verified for eventual problems. The report returned 5 problems, regarding the river connectivity. This is the exact message: End of checking river: oid=5 - PROBLEM: connectivity missing. Though I have tried solving this by merging the rivers (as I have merged some sub basins to better match the real situation of the watershed) or by manually snapping the rivers(it wasn't the case because they were already snapped to the main river, but just in case) and still I receive the same report with these problems. The hydro network is counting 12 rivers and 5 of them have this connection problem.

Has anyone faced this issue and solved it?


texts: Easy Extraction of features from a Feature dataset ( Not feature class)

texts:
My Feature dataset (residing inside a GDB) contains more than 30 feature classes and every month I like to extract features from this GDB based on one attribute, which is repeating in all the feature classes. For example Job_Month = "May".

Is there any tool available in ArcGIS which will extract all the features in a Geodatabase (or Feature dataset) at one time using an SQL statement? The extract data wizard does not have any filter/SQL for attributes.

Expected simple solutions, but seems not an easy one.


texts: Lauching pgsql2shp request on server by user

texts:
I am running PostgreSQL 9.4 and PostGIS 2.3.

I use Python scripting to launch a connexion to PostgreSQL with psycopg2.

My script is about the export of some queries in shapefile. Everything works fine when I launch it from my Windows session (I am administrator on the server), but when I try to launch it from another session I get an "Access denied" message. I am guessing it is linked with the fact that this other session isn't administrator.

I don't know much about systems, so I am guessing more informations are needed to help me solve that, but would anyone have an idea on where lies the problem and how to solve it? Obviously defining all my users as administrators is not an option.

Here is a sample of how my scripts look like (I connect through SSPI so don't need any user/password):

os.environ['PATH'] += r';SET PGCLIENTENCODING=LATIN1;S:\\ProgCommun'
os.environ['PGHOST'] = '10.159.107.13'
os.environ['PGPORT'] = '5432'
os.environ['PGDATABASE'] = 'SIEA'

cmds = '''"\\\\10.159.107.12\\xxx\\pgsql2shp.exe" -f "{}" -g geom DABASE "select * from table" '''
cmds = cmds.format(fichier)
subprocess.call(cmds, shell=True)

The IP are different since I try to use a local version of the pgsql2shp (and in so doing try to avoid executing anything on the server). But it doesn't work either way...


texts: Value labels not showing for raster when styled with PyQGIS

texts:
I'm developing a script which includes a function to style raster layers (based on  Setting 'classificationMode = Continuous' to a raster layer using PyQGIS? and Qgis: style a discrete raster). The script executes without any errors and colours the layer correctly. But the styled raster does not show any value labels (see image below).

If I go to the Style tab in Layer Properties and click apply, the values are then shown, so I think I must just be missing some sort of refresh command. My understanding was that lyr.triggerRepaint() should take care of this issue. I also tried iface.mapCanvas().refresh() but no such luck.

def styleRaster(lyr, lst):
    fcn = QgsColorRampShader()
    fcn.setColorRampType(QgsColorRampShader.INTERPOLATED)
    fcn.setColorRampItemList(lst)
    shader = QgsRasterShader()
    shader.setRasterShaderFunction(fcn)
    renderer = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1, shader)
    lyr.setRenderer(renderer)
    lyr.triggerRepaint()
    return

tciLyr = iface.addRasterLayer(hydroLyrs['tci'])

lst = [QgsColorRampShader.ColorRampItem(0, QColor("#f7fbff")),\
    QgsColorRampShader.ColorRampItem(4, QColor("#deebf7")),\
    QgsColorRampShader.ColorRampItem(8, QColor("#c6dbef")),\
    QgsColorRampShader.ColorRampItem(12, QColor("#9ecae1")),\
    QgsColorRampShader.ColorRampItem(16, QColor("#6baed6")),\
    QgsColorRampShader.ColorRampItem(20, QColor("#4292c6")),\
    QgsColorRampShader.ColorRampItem(24, QColor("#2171b5")),\
    QgsColorRampShader.ColorRampItem(28, QColor("#08519c")),\
    QgsColorRampShader.ColorRampItem(32, QColor("#08306b"))]

styleRaster(tciLyr, lst)


texts: How to use Local Azimuthal projection to convert Local NED to Geodetic LLA?

texts:
I'm trying to use Proj4j.EquidistantAzimuthalProjection to get the latitude/longitude of a new location based on its relative distance and bearing to a known location. Here is my source code:

import org.osgeo.proj4j.ProjCoordinate
import org.osgeo.proj4j.proj.EquidistantAzimuthalProjection

EquidistantAzimuthalProjection proj = new EquidistantAzimuthalProjection(43.694195, -79.262262)
proj.setUnits(org.osgeo.proj4j.units.Units.METRES)
ProjCoordinate src = new ProjCoordinate(10, 10)
ProjCoordinate dst = new ProjCoordinate
proj.inverseProject(src, dst)

Since I already set units to Meters, and 10m is a very small distance. The result should not be too far from latitude/longitude of the starting point  (43.694195, -79.262262). Hoevever the result yield by inverseProject is:

(138.66016814892942,-16.49593292366875)

The document clearly states that output unit of inverseProject is degree and input unit is consistent with the Projection itself (which is already set to meters). Do I miss something in invoking this library?


texts:
I'm creating a QGIS Python plugin and i want to execute some code after hitting the close button 'X'.
I've added the code bellow to init function and the closeEvent function but still not working

self.dlg.closeEvent = self.closeEvent


texts: How to access shapefile m-values with pyqgis?

texts:
How can I access (read/write) m values of a line Shapefile using pyqgis?

I know how to do it with arcpy cursors (http://resources.arcgis.com/en/help/main/10.1/index.html#/Reading_geometries/002z0000001t000000/ ) but I don't know how to do it with pyqgis.

Is it possible with the QGIS python library?


texts: Determine if a great circle is on the left or right from a point defined by latitude, longitude and bearing

texts:
Given

A point coordinate on the map specified by latitude and longitude. Shown as "C" on pictures below.
Direction where this point is moving, specified by bearing in degrees (0-359), where 0° means the point is moving strictly to North, 90° - to East, etc.
Two points (A, B) specified by latitude and longitude. These points define a line that point C should move towards.

Question
How can I determine whether the line is on the left or right side from the point C?
Examples

What may help
I think this problem can be solved using cross product calculation like in answers to this question, but I'm not sure how I would apply it having the coordinates defined by latitude and longitude.
I need to solve this problem using Swift programming language, but a solution using any other language or formulas would suffice.


texts: PostGIS function st_makevalid crashes the PostgreSQL Server

texts:
When I run my code running on PostGIS version 2.5, I don't get an error but when I run it on 3.1 it the server gives an error and enters the recovery mode every time.
pgadmin4 messages:(*"Connection to the server has been lost."*).

postgresql log : "*WARNING:  terminating connection because of crash of another server process
DETAIL:  The postmaster has commanded this server process to roll back the current transaction and exit, because another server process exited abnormally and possibly corrupted shared memory.
HINT:  In a moment you should be able to reconnect to the database and repeat your command.
WARNING:  terminating connection because of crash of another server process
DETAIL:  The postmaster has commanded this server process to roll back the current transaction and exit, because another server process exited abnormally and possibly corrupted shared memory.
HINT:  In a moment you should be able to reconnect to the database and repeat your command.
WARNING:  terminating connection because of crash of another server process
DETAIL:  The postmaster has commanded this server process to roll back the current transaction and exit, because another server process exited abnormally and possibly corrupted shared memory.
HINT:  In a moment you should be able to reconnect to the database and repeat your command.
CONTEXT:  SQL statement "SELECT (*

I have installed PostgreSQL on a separate server. Here I also get an error that PostgreSQL crashes due to the PostGIS function (eg:st_makevalid). In my opinion, there is a problem arising from the PostGIS version.
Can you try the queries yourself by trying on PostGIS extension versions 3.1.0 to 2.5?
I think we can see if there is a general problem of PostGIS 3.1.0 version.
PostGIS 3.1.0 - Query that raised an error:
select st_makevalid(
        st_force2d(
            st_multi(
                st_buffer(st_setsrid(ST_GeomFromGeoJSON(geometry), 4326), 0)
            )
        )::geometry (multipolygon, 4326)
    )
from (
        select j->'geometry' as geometry
        from (
                select (
                        '{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[39.75359,37.23446],[39.75359,37.23449],[39.75359,37.23415],[39.75391,37.23409],[39.7539,37.2341],[39.75373,37.23413],[39.75359,37.23446]]]}}'::json
                    ) as j
            ) t1
    ) t1

PostGIS 3.1.0 - Query that works without any problem:
select st_makevalid(
        st_force2d(
            st_multi(
                st_buffer(st_setsrid(ST_GeomFromGeoJSON(geometry), 4326), 0)
            )
        )::geometry (multipolygon, 4326)
    )
from (
        select j->'geometry' as geometry
        from (
                select (
                        '{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[28.89686,41.00149],[28.89698,41.00147],[28.89704,41.00157],[28.89693,41.0016],[28.89686,41.00149]]]}}'::json
                    ) as j
            ) t1
    ) t1

Operating System: CentOS Linux 8
CPE OS Name: cpe:/o:centos:centos:8
Kernel: Linux 4.18.0-193.28.1.el8_2.x86_64
Architecture: x86-64
postgresql-13
extension : postgis vers. 3.1.0


texts: Zoom to a collection of graphics in ArcGIS JavaScript API 4.x

texts:
Im trying to zoom to the extent of a couple of graphics that are on the map. None of the methods that worked for api 3.x are available in v4. Im hoping someone knows alternatives.

  https://developers.arcgis.com/javascript/latest/guide/functionality-matrix/index.html

graphicsUtils is marked as coming soon. This has the graphicsExtent(graphics) method that i could use.

I have tried the method

view.goTo()

but I need to specify the zoom myself. Which is not helpful. I want it to dynamically calculate bounds/extent and zoom to fit the graphics.


texts: Copy attributes from features in one layer in QGIS to another layer

texts:
I am building lake bottom topology maps for Lowrance HDS depth finders. Each iteration of map building (addition of new data) generates a new .shp file. To utilize the Lowrance tool that merges .shp files to a map (.at5 file), the features must have certain attributes defined and a value placed in the attribute.

Question is, once I have a .shp file with the attribute fields on the features, and a new .shp file is generated, how can I easily add the attributes to the new .shp features? Is there a way to copy/paste attributes? Or must I re-add each attribute to the new .shp file?

I have seen the postings about saving a feature in the old file and merging it with the new file, but it seems like there should be a simpler way.

Perhaps not. I have not been able to find a way.


texts:
I have three questions about spatial indexes management in GeoDjango, and I can't find the answers in the documentation.
First, I know GeoDjango creates a spatial_index for a geometry field by default when a table is loaded. Is the spatial_index automatically re-created when the table is updated? For instance, when new features are added to the table through the web application.
Second, If a table has more than one geometry fields, is GeoDjango creating spatial index for each of them or just for the first it encounters?
Third, I would like to create many spatial indexes on the same geometry table based on an id. Is it possible to do it with GeoDjango? For instance, many countries limits are stored in a single table and I would like to create a different spatial index for each country in the case each country are always queried individually.


texts:
I have a polygon layer and a point layer. According to the QGIS documentation, the childlayer has to have a foreign key. Unfortunately, my child layer has none. Is there a way to form the relation by geography only? So that even without having a proper foreign key, clicking on a polygon with the identify tool shows me all the points that are within this selected polygon in the attribute widget as child features?

I am aware I could create a new column in my child layer and assign that foreign key with the processing tools QGIS provides. As I have many layers to connect, this would be a rather time consuming process to do for all relations I want to establish!


texts: Accessing the output of a Python script tool dialog ArcGIS 10.1

texts:
I have a a fairly simple script tool dialogue which has been setup to get 6 parameters from the user and return them as the tool output as a list. the aims is that the out put would then be used elsewhere within a python add-in tool. I know I could revert to storing these values as a text file but that's not my preferred option. I can see the result in the results window but just can't get to it.

Below is the script to dialogue calls on closure and the seventh parameter is set to an output to store the first six parameters. I have been using the pythonaddins module to call the tool through PToolDialog(Toolbox,Toolname)

What I really want is to parse the output parameter from the script tool to the Python addin.

Script called by Python tool dialogue:

import arcpy

para1 = arcpy.GetParameter(0)
para2 = arcpy.GetParameter(1)
para3 = arcpy.GetParameter(2)
para4 = arcpy.GetParameter(3)
para5 = arcpy.GetParameter(4)
para6 = arcpy.GetParameter(5)
paraList = [para1,para2,para3,para4,para5,para6]
arcpy.SetParameter(6,paraList)


texts: Passing token to ArcGIS Server REST service?

texts:
I'm using .NET code to submit GET and POST requests to secured REST services served up by our internal ArcGIS Server 10.4.

How should I pass the security token?

Perhaps I can just add it as a querystring parameter but I do not know the parameter key.

I found documentation explaining the what and why of ArcGIS Server tokens but not how to actually pass them in a request. I have the url for obtaining a token. I just don't know how to use it.


texts: ArcGIS javascript Cluster Layer not working on a map with UTM spatial reference

texts:
I use this sample's clusterlayer module to create a cluster layer. problem is that when map spatial reference is UTM the cluster layer that I added to map does not work properly. the errors are
svg.js:43 Error: <circle> attribute cy: Expected length, "NaN".
setShape @ svg.js:43
svg.js:49 Error: <text> attribute x: Expected length, "NaN".
setShape @ svg.js:49
svg.js:49 Error: <text> attribute y: Expected length, "NaN".

But when the map's spatial reference is web mercator it works fine.
this is how i create my cluster layer:

var clusterLayer = new ClusterLayer({
  "spatialReference": map.spatialReference,
  "data": clusterData,
  "distance": 75,
  "id": "clusters",
  "labelColor": "#000000",
  "showSingles": true,
  "resolution": map.extent.getWidth() / map.width,
  "singleColor": pColor
});

I should note that clusterData is created based on map's spatial reference (coordinates are in map's spatial reference)

this is how clusterData created in UTM mode:

                        var splitedWkt = pLayerData[i]["SHAPE_WKT"].split(' ');
                        var endPoint = {
                            "x": splitedWkt[3],
                            "y": splitedWkt[4]
                        };
                        clusterData.push({
                            "x": endPoint.x,
                            "y": endPoint.y,
                            "attributes": {}
                        });

"SHAPE_WKT" is the feature's shape that i get directly from database. (Well Known Text).

And this is one of the objects in clusterData array:

{x:"565623.31040000",
 y:"3616457.34120000",
 attributes:{}

}


texts: Why does ArcGIS Desktop 10 HTML popups give JavaScript/HTML Error?

texts:
Using ArcMap 10 Desktop HTML Popups.

The issue I am having is that the html < select > elements aren't working - you can select it and use the arrow keys to navigate through(up/down) the list but you can't bring up a drop down with the visible list bit.

When on desktop IE running in ie 7 or 8 mode selects don't work I get a popup notification but it wont allow me to unblock it.

This a bug, is there a fix? can I upgrade the IE higher than 7/8? chromeframe? I just need to get the selects working.

< ? php

#$tmp = "< select id='something' >";

for(something){

$value = "value";

$name = "name";

#$tmp.="< option value="$value" >$name< / option >";

    }

#$tmp.="< /select>";

echo $tmp;

? >


texts: Geodesic distance discrepancies using proj

texts:
I have two positions represented as latitude,longitude pairs, and I'm trying to get the distance between them. I've 3 methods that all yield different results, where I think they should be VERY close. The methods are:

Converting the lat/lon to a Cartesian UTM system, and then using the Cartesian distance
pyproj.Geod
The "geod" command-line tool

For an arbitrary pair of points that are ~1km apart, I see:

diff utm:          1161.43346448816
diff geod:         1161.1017567084384
diff geod cmdline: 1163.971

This looks wrong. The UTM probably should be slightly off because the two Cartesian systems are slightly tilted, but is 30cm over 1km really how far off it should be?

The difference between pyproj.Geod and the "geod" cmdline tool makes even less sense. In theory they're doing the same thing, so the results should be identical.

A Python script to run the 3 methods, and print the ranges:

#!/usr/bin/python3

import sys
import pyproj
import numpy as np
import subprocess

lat0,lon0 = (31.4617, -83.4849)
lat1,lon1 = (31.4704, -83.4781)

zone = 1 + int((lon0 + 180) // 6)
print(f"zone: {zone}")

toutm = \
    pyproj.Proj(proj           = 'utm',
                zone           = zone,
                ellps          = 'WGS84',
                preserve_units = False # output in meters
                )

e0,n0 = toutm(lon0, lat0)
e1,n1 = toutm(lon1, lat1)

print(f"diff utm:          {np.sqrt((e0-e1)**2. + (n0-n1)**2.)}")
print(f"diff geod:         {pyproj.Geod(ellps = 'WGS84').inv( lon0,lat0,  lon1,lat1)[2]}")

cmd_input = f"{lat0} {lon0} {lat1} {lon1}"
cmd_output_range = subprocess.check_output( ("geod", "+ellps=WGS84", "-I", "+units=m"),
                                            input = cmd_input,
                                           encoding = "ascii").split()[-1]
print(f"diff geod cmdline: {cmd_output_range}")


texts:
I have a number of aerial images that do not align with where my GPS says I am. For example, I am standing at the intersection of two roads with a Trimble R10  and a handheld IOS device with the aerial imagery as a geotiff displaying on it. The highly accurate GPS position is saying I am in the ditch about 6m away on the imagery. To rectify the imagery alignment issue I seem to be having I would imagine some program of placing Ground Control Points is necessary. My question is, how do you do this for aerial imagery that is already existing, not being flown at time of GCP collection. I know I need to find easily identifiable features on the aerial image but (1) I am in the middle of no-where and easily identifiable features are at a premium and (2) I will have to manually choose where I think the feature is on the image which could be off by up to, even over a meter thus defeating my Trimble's high accuracy. I think I am approaching this the only way you can, but want to see if there are different thoughts on the topic.

Trimble R10 Model 2 - (https://www.duncan-parnell.com/storage/attachment/datasheet_-_trimble_r10_model_2_gnss_system.pdf)


texts:
I have set up a page layout using the Altas Feature in QGIS v3.2.3. The layout has been configured to show the location of a polygon and the location of vertices that make up each of these polygons. Alongside this I have added an attribute table to list the X and Y values of each vertex for every polygon. The Atlas feature itself is linked to the polygon shape layer whilst the attribute table is linked to a vertex point layer. And whilst there have been no problems exporting this set-up there is one aspect which I cannot seem to resolve which is that in some instances the attribute table lists vertices from different polygons (rather than listing vertices only for the Altas Feature). I have tried toggling with some of the item properties for attribute table e.g. selecting 'show only features intersecting altas feature' however I believe that the solution lies with the expression filter. Unfortunately I have not been able to set up an expression filter that will work with the Atlas Feature and therefore I am looking for some advice on this.

I have attached an image for reference to the above.


texts: How to use optional inputs in model processor for QGIS 2.18/3.0?

texts:
I have developed a processing model in QGIS 2.18 where all my inputs are optional and I need them to be truly optional. However, the model will not run unless I have all inputs present.

I have tried this model in both 2.18 and 3.0. I would prefer to use 3.0, but cannot figure out how to edit the python script.  I think my problem is that it is not registering outputs as optional, but I'm not sure. Below is a sample of my code from the 2.18 model. I have also attached a screen capture of the start of the model in 3.0. I have been successful at running the models if all inputs are entered. I've tried recreating the model on two different computers, tried using QGIS 2.14 and searched the stackexchange/google etc.

##lines to points=name
##accesspoint=optional vector
##cableburied=optional vector
##ugspan=optional vector
##equipment=optional vector
##attachment=optional vector
##buildingfootprint=optional vector
##spliceclosure=optional vector
##buriedugspan=output vector
##ugspanvaccesspoint=output vector
##cableburiedv=output vector
##splicecable=output vector
##ugspanv=output vector
##ugspanaccess=output vector
##buildingcable=output vector
outputs_SAGACONVERTPOLYGONLINEVERTICESTOPOINTS_1=processing.runalg
('saga:convertpolygonlineverticestopoints', cableburied,cableburiedv)
outputs_QGISJOINATTRIBUTESBYLOCATION_3=processing.runalg
('qgis:joinattributesbylocation',
spliceclosure,outputs_SAGACONVERTPOLYGONLINEVERTICESTOPOINTS_1['POINTS'],
['overlaps'],0.0,1,'sum,mean,min,max,median',0,splicecable)
outputs_SAGACONVERTPOLYGONLINEVERTICESTOPOINTS_2=processing.runalg
('saga:convertpolygonlineverticestopoints', ugspan,ugspanv)
outputs_QGISJOINATTRIBUTESBYLOCATION_4=processing.runalg
('qgis:joinattributesbylocation',
equipment,outputs_SAGACONVERTPOLYGONLINEVERTICESTOPOINTS_1['POINTS'],
['overlaps'],0.0,1,'sum,mean,min,max,median',0,None)
outputs_QGISJOINATTRIBUTESBYLOCATION_5=processing.runalg


texts: Building exaggerated elevation profiles from VRT DEM in QGIS

texts:
I have this little challenge, where I need to build an elevation profile from a VRT DEM mosaic of +5000 TIF files. All Z-values should be exaggerated by factor 10 as it's common in road construction.
The profile must be exported to DXF for further processing in AutoCAD.
Currently I'm using the "Profile Tool" plugin to produce the profiles, but there is no option to exaggerate the output.
What's the straight forward method to complete this task in QGIS?


texts: Have GraphicsLayer working but now need to put it in or convert to FeatureLayer (JavaScript)

texts:
I have a GraphicsLayer working for "Regions" and the corresponding popup template information for the staff within each region (below).

        function fillRegionMarkers() {

            // Create a symbol for drawing the point
            var markerSymbol = {
              type: "picture-marker",  // autocasts as new SimpleMarkerSymbol()
              url: "https://SomeUrl/map-pin-16.png",
              height: "14pt",
              width: "14pt"
            };

            // Create Graphic Layer
            var gl = new GraphicsLayer({ id: "RegionPoints" });

            //outside loop: REGIONS
            $.each(arrRegions, function(idxRegion, elementRegion) {

                // First create a point geometry
                var point = {
                  type: "point",  // autocasts as new Point()
                  longitude: elementRegion.Longitude,
                  latitude: elementRegion.Latitude
                };

                var DCODCEInfo = "";

                //inside loop: DCO/DCE STAFF
                $.each(arrDCODCE, function(idxDCODCE, elementDCODCE) {

                    if (elementRegion.Id == elementDCODCE.Region.Id) {
                        //Should never be null, so no need for... if (elementData.lat && elementData.lon)
                        //lat and long should never be both 0, so no need for... if (elementData.lat != 0 && elementData.lon != 0)
                        DCODCEInfo += "Position: " + elementDCODCE.Position + "
" + "POC: " + elementDCODCE.StaffInfo + "

";
                    }
                })

                var template = {
                    title: "
(All Regions Map)

Region: " + elementRegion.FEMARegion + "
Location: " + elementRegion.Location,
                    content: [{
                        type: "text",
                        text: DCODCEInfo
                    }]
                };

                // Create a graphic and add the geometry and symbol to it
                var pointGraphic = new Graphic({
                  geometry: point,
                  symbol: markerSymbol,
                  popupTemplate: template
                });

                //layer.add(pointGraphic);
                gl.graphics.add(pointGraphic);
            })

            //Add graphics layer to map
            map.add(gl);
        }

The problem is a URL feed with a FeatureLayer is now being used to show Regions (below).

    var featureLayer = new FeatureLayer({
      url:
        "https://SomeUrl/Regions/FeatureServer/0",
      title: "Regions",
      visible: true

    });

    map.add(featureLayer);

This FeatureLayer does not have the DCO/DCE staff or pinned lat/long locations like the GraphicLayer has. So how do I add or convert this GraphicLayer to a FeatureLayer?


texts: Order a list of points to make a polygon with smooth outer border

texts:
I have a function that pulls some lat/longs from a Database and then I put those values in a polygon to see it's shape. The problem is they are so random in position and direction from each other that I get a ton of triangle shapes due to the overlapping. Is there a way to order the positions so I get a solid polygon without the overlapping holes?
 Tried:

Dim orderedPoints = points.OrderBy(Function(ll) ll.Lng).ThenBy(Function(ll) ll.Lat).ToList

I feel it needs to be drawn in a circular pattern, finding the farthest as it moves around in a given direction. Hoping to not have to reinvent the wheel if poss here. If there is a name for this process I would be more than happy to research it once I know what it is called.

  EDIT: Found this for Convex Hull. Made this code so far. Not sure how to proceed from here.

Public Class ConvexHull
   Public Function ccw(a As PointLatLng, b As PointLatLng, c As PointLatLng) As Integer
     Dim area2 = (b.Lat - a.Lat) * (c.Lng - a.Lng) - (b.Lng - a.Lng) * (c.Lat - a.Lat)
     If (area2 < 0) Then Return -1
     If (area2 > 0) Then Return +1
     Return area2 ' basically collinear
   End Function
   Public Function collinear(a As PointLatLng, b As PointLatLng, c As PointLatLng) As Boolean
     Return ccw(a, b, c) = 0
   End Function
 End Class


texts: Removing RubberBand item in QGIS 2.18.19

texts:
I have QGIS 2.18.18 and 2.18.19 in seperate computers. I use the following script to remove RubberBand items. It works in 2.18.18 but doesn't work in 2.18.19

import qgis

# Get RubberBands
rbs = [ i for i in iface.mapCanvas().scene().items()
            if issubclass(type(i), qgis._gui.QgsRubberBand) ]

# Remove RubberBands
for rb in rbs:
    if rb in iface.mapCanvas().scene().items():
        iface.mapCanvas().scene().removeItem(rb)

Colored polygons are RubberBand items.


texts: Getting latest image from imageCollection of Google Earth Engine

texts:
I'm trying to get the last image available in an image collection but I can't find any function that works.

I've tried also to sort the collection in the reverse order, but it doesn't work.

here is my code:

var chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/PENTAD");

//get the first image and its date
var firstimg = chirps.limit(1, 'system:time_start').first();
var firstdate = ee.Date(firstimg.get('system:time_start'))
print(firstdate);

//now I try to get the latest...

chirps.sort('system:time_start', false);
var lastimg=chirps.limit(1, 'system:time_start').first();
print(ee.Date(lastimg.get('system:time_start')));

How can I get the last image available and its date in order to know which is the last date available?


texts: Copy>Paste from meter to feet coord-sys, and coordinates don't change in ArcMap 10

texts:
My dataframe coordinate system is in UTM with linear unit Meter.

If I have a Polygon shapefile in a UTM coordinate system with a linear unit of Foot_US and I am sent a shapefile in the same UTM coordinate system with a linear unit of Meter.  I reluctantly copy/paste the linear unit Meter feature into the linear unit Foot_US shapefile and the vertices' xy are identical in both the linear unit Foot_US and Meter.

What am I missing.  I would have figured they would be different.


texts: ArcGIS 9 doesn't recognize ArcGIS 10 .lyr files

texts: I symbolized a parcel shapefile in ArcGIS 10 with about 2 dozen different colors. I need to use this in ArcGIS 9.3 as well, so I saved it as a .lyr file. This works fine in 10. But when I go to 9.3, I can add the shapefile, but not the .lyr file, or import it for symbolization (geometry doesn't match error). Is there a way to make a .lyr file in 10 so that 9 will recognize it? The .lyr file comes from the same shapefile, so I don't know why it wouldn't recognize it. I want to use this symbolization for a half dozen shapefiles a/o feature classes, so really don't want to manually recreate it. Thanks.


texts:
I am installing Whitebox GIS software on linux.

I have already:

set up a clean xfce4 machine
installed java 8 successfully
unzipped the latest version of whitebox (Iguazu 3.2.1) to the desktop

Whenever I click on the JAR file or try to run through  terminal nothing happens.

Am I missing some key dependencies?


texts: How to compare schemas from two File geodatabases?

texts:
We have template geodatabase, which we send out to our regional Data creators. They update these geodatabases and digitize all the data from their region in the predefined featureclasses.

Now that we have got the File geodatabases, we find that often many of the featureclasses have been modified. Maybe a field might be added or removed. In some cases, the Feature-class itself is deleted or a new one created.

I would like to get a report indicating which feature-classes and tables are changed, as well as the differing Fields in the common tables & featureclasses.

I have already looked at several questions such as:

Tools to compare unversioned copies of geospatial data
Is there a good way to get a diff between an old and a new version of an SDE database?

But the answers given in these Questions have not been helpful.


texts: Problem with projection in meters with ArcGIS 10.1

texts:
I hope you can help with the following problem.
I have a shapefile and want to project it using a projected coordinate system with meters as unit. What I did was to change the projection from ArcCatalog and then importing the shapefile into ArcMap. But when I do this, there are two problems:

Distances are completely wrong. By using the Measure tool, houndreds of kilometers are measured as less than a meter!
The numbers on the bottom-right which indicated the position and the units of the map are also wrong. When displayed in Decimal Degrees, are either 0 or 1. When I change this in meters, numbers are still very small.

The problem persists even after the projection is added to the dataframe. What am I doing wrong? Thanks for any help and suggestion.


texts: GridCoverage2DReader not reading 16 bit grayscale?

texts:
I have a geotiff that has 16-bit gray scale pixels and GridCoverage2DReader instance only returns the nodata value. Yet both QGIS and uDig show the image and I see WGS 84 and EPSG:4326 for the geotiff properties.

I am using geoserver 2.8.2 with the DDS/Bil extension and Java 1.8.x.

Any ideas what might be wrong?

I already have debug enabled for logging and did not see any exceptions thrown but did not look for debug logging for the reader. I will look again for relevant logging. I am trying to find the implementation class for the reader since I suspect what I mentioned is probably just an interface. Right now all I know is the reader is coming from a resource pool. I think I did notice one of the classes would let me get to read parameters for the reader. I will get that info and see if it makes sense. I suspect that the reader was setup by GeoServer and not the bil plugin since it has the store info and layer info for the GeoTIFF file.

I am a newbie to GeoServer but not to java.


texts: Only displaying polygons over certain size in ArcMap using Definition Query?

texts:
I have a polygon containing lots of bodies of water and only want to display the bodies of water over a certain size.

Can I do this using a definition query?


texts: Clip Point Features to Polygons with ModelBuilder/ArcPy

texts:
I have point features for about 2000 individual trees spread throughout 50 different tree districts. The tree districts are polygons in a polygon feature class and the trees are point features in a point feature class.

I would like to create a mapbook showcasing each individual tree district and the corresponding trees in each district. I could use Data Driven Pages, but I need closeups of each tree district (polygon) removed from its surrounding districts.

I could accomplish this by running a definition query to display a single tree district, clip the trees to that district, and then create a mapbook page. But I'd have to do this 50 times.

How can I use ModelBuilder or ArcPy to perform the aforementioned process?

I've tried using How to Batch Clip in ArcGIS Desktop 10 using Python/ArcPy? as a guide but it doesn't exactly fit my needs.


texts: Running weekly sceduled project with different input data on QGIS Server

texts:
I would like to run a QGIS project with a QGIS server, weekly scheduled. I wish the project to run with a different input data, a shapefile layer.
Could such a setup be automatized with QGIS server?


texts: Leaflet Vector Layers Not Rendering During FlyTo Animation

texts:
In version 1.0 of Leaflet, when using the FlyTo function as per this JSFiddle, is it possible to force vectors to render/ draw as the animation is happening?

My current issue is when using the marker.cluster function, the animation starts and then none of the points display as you would expect until the animation ends.

I believe the issue is like this OpenLayers question, is there a similar solution for Leaflet?

UPDATE: It appears the issue is related predominately to GeoJSON layers, see this first example that animates well and then this example that uses a GeoJSON layer and the animation doesn't render the vector.


texts: ArcGIS Javascript Editor - showAttributesOnClick feature

texts:
I'm using the Editor dijit to add feature to my map. I was looking at one of the online samples, which mentions a "showAttributesOnClick" setting in the description. When I look at the code in that same sample, that setting isn't used and it's not documented in the API. Does that feature still exist? Has it been replaced with something else?

I'd like to only display the AttributeInspector for a new feature, and hide it for existing ones.

Thanks for the help!


texts: Filegeodatabase fanout and Python problem

texts: I have fme wb which does fanout with gdb. It means that in the end as output sets of data I get few new created gdb. But If I start this fme wb within python script It does not do fanout. Just write data in one gdb. i do not understand where could be a problem?


texts: trouble with i.landsat.toar failing to identify satellite

texts: I am trying to use the i.landsat.toar tool on L7 bands. I have tried to use both the MTL and MTLold, but no matter what I try, it returns "Failed to identify satellite." What am I doing wrong?


texts:
I know that there are many scripts (mostly in matlab) out there that address the spread of disease/zombie infection, but does anyone know of ways this has been developed/implemented in ArcGIS?

Essentially, how would you deal with taking into account issues of:

Terrain
Population density (by polygon or raster)

Would you write it in something like python or are there tools/models/scripts that already exist that I don't know about?

From what I can find, a majority of the modelling is done in math, like here.

What steps would you take personally to answer the question: Given point X, how fast would a zombie outbreak spread given Y speed, Z terrain, or other variables?

The reason I ask, I teach GIS to a lot of young people and would love to use this as an opportunity to inspire them.

Thanks!


texts: Leaflet filter on geoJson properties : multiple criteria

texts:
I'm using Leaflet filter on geoJson data. I'm able to apply it with a simple criteria :
filter: function(feature, layer) {
     return feature.properties.field1 == "yes";
}

I can also apply two criteria with logical operators like && or || :
filter: function(feature, layer) {
     return (feature.properties.field1 == "yes" && feature.properties.field2 == "yes");
}

I want now to insert a condition (depending of a variable):
filter: function(feature, layer) {
   if (varState == "01") {
     return (feature.properties.field1 == "yes" && feature.properties.field2 == "yes");
   }
   else if (varState == "02") {
     return (feature.properties.field1 == "yes" && feature.properties.field2 == "no");
    }
    else {
      return feature.properties.field1 == "yes";
    }
}

Does someone has already succeed this ?


texts:
I have 5 historical aerial sets for my county each consisting of several hundred SID files. Previous employees have created Image Catalogs for each. Now I want to turn these into services I can share with ArcGIS for Server 10.1. Running the mosaic tool is taking forever, I started at 830am and as of 4pm it is only through 30 of 508 files. There has to be an easier way?

Is there a workflow or a way that will be quicker?


texts: Most concise way to check existence of ALL files in a list (ArcPy)?

texts:
I'm running a process that includes generating a set of 3 curvature layers. For various reasons, sometimes I may generate just 1, 2, or all 3. Later I may need to go back and generate another layer, but the complication is that to create the first set of layers I used a randomly generated temporary surface. So if I were to create a sibling curvature layer, it wouldn't match the original. What I do, then, is to see if any of the 3 layers are missing. If so, I delete and recreate them all 3. But there must be a cleaner way.

curvLayers = [curvpath, profpath, planpath]
if (not arcpy.Exists(curvpath) or
  not arcpy.Exists(profpath) or
  not arcpy.Exists(planpath)):
    for layer in curvLayers:
        if arcpy.Exists(layer): DM.Delete(layer)
        <create new layers>

I'd love to think there was something like:

if not all arcpy.Exists(curvLayers):


texts: QGIS increase number of recent selections

texts:
is there any possibilty to increase the number of recent selections in the menu select by expression?
I can´t store more than 20 selections. If I create a new selection, the oldest selection is overwritten.

Or is there a way to create a own group menu, like Recent(Selection), with a own description and selection?
I use the QGIS Version 3.4.8 LTR


texts: Joining features based on location and attribute at the same time

texts:
I have a task of assigning data from database that I represent as point layer, and containing information about characteristics of soils through the world, and also soil type. I also have polygon layers for each lattitude, that contain area of each type of soil (same ones as point layer, therefore can be ussed as join feature). Inside each lattitude layer, each soil type is only one, multipart polygon. I need to assign values from point layer in a way that the soil type always match, but also that it uses point closest geograpicly. If I use just join by attributes, it doesnt take into account spatial location, and vice versa, and I really have trouble figuring how I could approach that. Each layaer has around 40 soil types, and so doing it manually would take ages. I can also add that majority of points are way off the polygons that I intend them to join to, but as long as the soil types match and points is closest, its what I need.
I have basicly no knowledge about python scripting.
 I tried making model to handle that, but its also beyond me.
Thanks in advance


texts: Custom data entry form (.ui) crashes qgis

texts:
I'm using Linux Mint 17.2, qgis 2.8.1 and have also tried this on 2.10:

Made a simple form using qt4 designer
The field names match those in the attribute table
used the properties page of qgis to select the "provide ui" option
click Apply and okay
using the attribute identify tool clicked on the layer

instantly crashes qgis.

BTW, this is my first question, but have learned a lot from those posted by others.


texts: OGR Python: Copying features from one datasource to another works with shapefile but not with geopackage

texts:
I have a shp/dbf EPSG4674 input datasource, and I want to copy features from one layer of this input datasource, create a new output datasource and write these features to that output datasource.

As of now, I am copying all features from one layer to the other, but later I will conditionally filter some of them and write to different datasources depending on the feature's attributes.

When my output datasource is a shapefile, it works perfectly:

from osgeo import ogr

ogr.UseExceptions()
input_datasource = ogr.GetDriverByName('ESRI Shapefile').Open(r'C:\Users\rodri\Documents\test\test_line.shp')
input_layer = input_datasource.GetLayerByIndex(0)

# UNCOMMENT
#output_datasource = ogr.GetDriverByName('GPKG').CreateDataSource('output.gpkg')

# COMMENT
output_datasource = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource('output')

output_layer = output_datasource.CreateLayer('output_layer', input_layer.GetSpatialRef(), input_layer.GetGeomType())

# Copying the old layer schema into the new layer
defn = input_layer.GetLayerDefn()
for i in range(defn.GetFieldCount()):
    output_layer.CreateField(defn.GetFieldDefn(i))

# Copying the features
for feat in input_layer:
    output_layer.CreateFeature(feat)

When I change the output datasource to a geopackage (comment the line below "# COMMENT" and uncomment the line below the "# UNCOMMENT" line), it no longer works and it gives me this error:

Traceback (most recent call last):
  File "gpkg_test.py", line 24, in <module>
    output_layer.CreateFeature(feat)
  File "C:\Users\rodri\Anaconda3\lib\site-packages\osgeo\ogr.py", line 1406, in CreateFeature
    return _ogr.Layer_CreateFeature(self, *args)
RuntimeError: failed to prepare SQL: INSERT INTO "output_layer" ( "fid", "", "ID", "AREA_ID") VALUES (?, ?, ?, ?)

(my layer is a multilinestring layer with two fields, ID and AREA_ID)

I think it's important to mention that my drivers are okay, I have downloaded an example script that generates a geopackage with both vector and raster data and it worked without any problems.

I'm running Python 3.6.5 on Windows with GDAL 2.3.2.


texts: Python code error when making fusion of multispectral image and DEM

texts:
I tried to fuse 4 multispectral bands from Landsat 8 with resampled SRTM, Aspect, Slope. It should contain 7 bands: R,G,B,NIR,DEM,Aspect,Slope. The code execution generates output raster, MS bands are ok, while DEM,Aspect and Slope are empty. Code results an error. What is wrong?

import rasterio
import numpy as np

file_list = ['E:/Burabay2019/SpaceImages/Calculated/LC08_L1TP_156023_20150626_20170407_01_T1_clip/LC08_L1TP_156023_20150626_20170407_01_T1_B4_c.TIF',
         'E:/Burabay2019/SpaceImages/Calculated/LC08_L1TP_156023_20150626_20170407_01_T1_clip/LC08_L1TP_156023_20150626_20170407_01_T1_B3_c.TIF',
         'E:/Burabay2019/SpaceImages/Calculated/LC08_L1TP_156023_20150626_20170407_01_T1_clip/LC08_L1TP_156023_20150626_20170407_01_T1_B2_c.TIF',
         'E:/Burabay2019/SpaceImages/Calculated/LC08_L1TP_156023_20150626_20170407_01_T1_clip/LC08_L1TP_156023_20150626_20170407_01_T1_B5_c.TIF',
         'E:/Burabay2019/SpaceImages/Calculated/DEM_clip/srtm_burabay_res30_c.TIF',
         'E:/Burabay2019/SpaceImages/Calculated/DEM_clip/srtm_burabay_res30_aspect_c.TIF',
         'E:/Burabay2019/SpaceImages/Calculated/DEM_clip/srtm_burabay_res30_slope_c.TIF']

with rasterio.open(file_list[1]) as src0:
   meta = src0.meta
meta.update(count = len(file_list))
with rasterio.open('E:/Burabay2019/SpaceImages/Calculated/datafusion_LC08_L1TP_156023_20150626_dem_20190728.TIF', 'w', **meta) as dst:
for id, layer in enumerate(file_list, start=1):
    with rasterio.open(layer) as src1:
        dst.write_band(id, src1.read(1))

The error is below:

ValueError                                Traceback (most recent call last)
<ipython-input-2-e20e01c8ca7c> in <module>
 24     for id, layer in enumerate(file_list, start=1):
 25         with rasterio.open(layer) as src1:
 ---> 26             dst.write_band(id, src1.read(1))
 rasterio\_io.pyx in rasterio._io.DatasetWriterBase.write_band()
 rasterio\_io.pyx in rasterio._io.DatasetWriterBase.write()
 ValueError: the array's dtype 'float32' does not match the file's dtype 'uint16'