texts: Mapbox GL JS set Max Zoom level not working?

texts:
I cant seem to get my map to set a max zoom level.

Here is a snippit of my code -

var map = new mapboxgl.Map({
container: 'map',
style: 'mapbox://styles/cdickinson11/cjg49l1f305ou2rls3hd65a6s',
center: [-77.034084142948, 38.909671288923],
zoom: 10
options.maxZoom: 13
});

I also tried this, adapted from here. It did not work either, I'm guessing because it is from mapbox js instead of mapbox gl js.

var options = {
maxZoom: 13
};
mapObject = L.mapbox.map('mapDiv', undefined, options);

I'm fairly new to this type of thing. I'm pretty sure my issue lies in a lack of basic understanding about javascript, and about the difference between mapbox gl js and mapbox js. Any help?


texts: QGIS3 in custom application via PyQGIS (Linux)

texts:
I am trying to establish a QGIS3 custom application under Linux (Arch Antergos), according to the documentation (section 1.4.2) the following should work (but it doesn't). I executed the subsequent code from the system command line:

export PYTHONPATH=/usr/share/qgis/python
export LD_LIBRARY_PATH=/usr/lib
python3
from qgis.core import *
QgsApplication.setPrefixPath("/usr/bin/qgis", True)
# running showSettings() also does not return the expected results
# QgsApplication.showSettings()
# running QgsApplication.setPrefixPath("/usr", True) looks better with showSettings(), however, fails also in the next step
qgs = QgsApplication([], True)
qgs.initQgis()

After running qgs = QgsApplication([], True) Python is terminated with following error message:

malloc_consolidate(): invalid chunk size
[1]    5581 abort (core dumped)  python3

The same code worked with QGIS 2.18.

Has something changed or am I missing something?


texts:
I have a SRTM DEM and I want to create a shaded relief from it. I created the shaded relief in GRASS and the result is very nice, but a little rough because the area is near flat and the DEM is 90m in resolution.

What I want is to make the DEM smoother in order to generate a smooth shaded relief. Is there an algorithm or interpolation method to do that?

Here is the shaded relief to get an idea, I want to flatten these small bumps:


texts: Calculate max distance area nodata cells in raster

texts:
I am aiming to write a tool in Arcpy that merges (mosaic) multiple raster files, fills the nodata cells (via gdal.FillNodata) and applies a hillshade to the raster. It is working fine but I want to optimize my script.

The raster files that are put in to the tool contain nodata cells that the tool fills with GDAL.FillNodata. However the sizes of area's that have these nodata 'gaps' vary a lot from raster to raster, so setting a maximum search distance (maxSearchDist) is quite a hassle. And for optimization purposes I don't want to set an extreme distance.

So my goal is to calculate the maximum cell distance of the biggest area with nodata cells within a set of raster files. And that'd be the input as my maxSearchDist. How would this be achieved? I am looking at raster calculations but I don't know whether to look for my answer in GDAL or Arcpy itself.

In the figure below I tried to visualize what I want to achieve. On the left is the raster file. From this raster I want to calculate which area of nodata cells is the largest and secondly (on the right indicated with the blue arrow) what the maximum distance of nodata cells is within that area.


texts: Moving vectors to specified coordinates in QGIS?

texts:
I have a layer of vectors located around, say, (100, 100), and I want to move them to new coordinate like (1000, 1000).

How should I do it using QGIS?


texts:
After importing PyQt4.QtCore,PyQt4.QtGui and processing and after having correctly defined layers, I executed this script in PyQGIS but it returns this error:
Traceback (most recent call last):
  File "", line 9, in
NameError: name 'height' is not defined (so it refers to line "print height").
I also failed by checking in previous posts:
like https://stackoverflow.com/questions/20177159/python-nameerror-name-width-is-not-defined or NameError: variable not defined even though set multiple times and other ones but do not fit because different from my case.
Any suggestions?

layerbld.startEditing()
maxh_field=layerbld.dataProvider().fieldNameIndex('MaxHeight')
for feat in layerbld.getFeatures():
  for f in layer.getFeatures():
     if feat['CR_EDF_UUI']==f['CR_EDF_UUI'] and feat['maxUN_VO_1']==f['Shape_Area']:
       height=f['UN_VOL_AV']
     else:
       pass
  print height
  layerbld.changeAttributeValue(feat.id(),maxh_field,height)
layerbld.updateFields()
layerbld.commitChanges()

I tried to correct in this way but it returns only null values, that's not possible:

for feat in layerbld.getFeatures():
  for f in layer.getFeatures():
    if feat['CR_EDF_UUI']==f['CR_EDF_UUI']:
      if feat['maxUN_VO_1']==f['Shape_Area']:
          height=f['UN_VOL_AV']
   else:
      height=0
  print height


texts: Seeking historic aerials for some areas in Dona Ana County, New Mexico?

texts:
How do I find several historic aerials from the mid-1970's or early 1980's and where can I download them?

I would like them to be free of charge.


texts: How can I copy non-blank values to a new column from different columns in the same layer in QGIS?

texts:
My data:

n1    n2    n3    NEW_COLUMN
11    null  null
null  33    null
null  null  55

My desired result:

n1    n2    n3    NEW_COLUMN
11    null  null  11
null  33    null  33
null  null  55    55


texts: Resampling Sentinel-2 ImageCollection on GEE having the bands with the same spatial resolution

texts:
I'm a beginner at Google Earth Engine.
I've seen on previous topics that resampling can be made on single images, like in this topic on downscaling raster images. I have an ImageCollection of Sentinel-2 data and would like to resample at least the bands with 20 m spatial resolution to a 10 m. This is the part of the code working on data collection.
var area = ee.Geometry.Polygon([[8.647458, 44.659079],[8.588429, 44.658882],[8.587969, 44.723788],[8.647064, 44.723985]]);

/**FUNCTION to MASK CLOUDS using the Sentinel-2 QA band
@param {ee.Image} image Sentinel-2 image
@return {ee.Image} cloud masked Sentinel-2 image**/
function maskS2clouds(image) {
    var qa = image.select('QA60');
    //Bits 10 and 11 are clouds and cirrus, respectively
    var cloudBitMask = 1 << 10;
    var cirrusBitMask = 1 << 11;
    //Both flags should be set to zero, indicating clear conditions
    var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));
    //Return masked and scaled data, without the QA bands
    return image.updateMask(mask).divide(10000).select("B.*").copyProperties(image, ["system:time_start"]);
}

//define collection parameters
var start = '2017-11-01';
var end = '2018-10-31';
var cld_max = 20;   //max % cloud cover allowed in image collection; for data until now used: 20

//collect S2 data
var s2 = ee.ImageCollection('COPERNICUS/S2_SR')
              .filterBounds(aoi)
              .filterDate(start,end)
              // Pre-filter to get less cloudy granules
              .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',cld_max))
              .map(maskS2clouds);

But, when it comes to the resampling procedure, doesn't seem to work
var b8 = s2.select('B8'); // 10-m resolution
var b8a = s2.select('B8A'); // 20-m resolution

// Get information about the S2 projection at different resolutions
var s2_10m = ee.Image(s2.first().select('B8'));
var proj = s2_10m.projection();
print('S2 projection (10m-res):', proj);

var s2_20m = ee.Image(s2.first().select('B8A'));
var proj = s2_20m.projection();
print('S2 projection (20m-res):', proj);

var transform = [10, proj.transform[1], proj.transform[2], proj.transform[3], -10, proj.transform[5]];
print('new CRS transform:', transform);

// Display resampled image collection
var b8a_res = b8a.first().resample('bilinear').reproject({crs: proj.crs, crsTransform: transform});
Map.addLayer(b8a_res, {'min': 0, 'max': 0.5}, 'nNIR_resampled');

And reports the

nNIR_resampled: Layer error: Can't encode object: function standin(_0) {
return newF.apply(this, arguments);
}

Is there a better way to resample the image collection 'B8A' (narrow NIR) to a 10-m spatial resolution (just like NIR and RGB bands)?


texts: Multivalue parameter as text in Python toolbox - reading in to imported module function

texts:
I'm building a Python toolbox with a tool which calls specific functions from a custom module.

The main input parameter for this tool is a multivalue list of raster layers.

In the execute block of my tool, I use parameter[2].valueAsText to read in said parameter, as normal, and then feed it into a custom module function as follows:

output          = parameters[0].valueAsText
input           = parameters[2].valueAsText
function        = parameters[1].valueAsText
functionname    = "Local" + function
option          = parameters[3].valueAsText
expression      = parameters[4].valueAsText

newGrid = getattr(localfuncs,functionname)(input,option,expression)
newGrid.save(output)

In this case, I've tested it with one custom function (corresponding to 'functionname'):

def LocalMean(input,option,expression):
    import arcpy
    newGrid = arcpy.sa.CellStatistics(input,"MEAN")
    return newGrid

The trouble is that the 'input' (multivalue read in as text) ends up being a semicolon-delimited string of file names, and the function in the custom module can't seem to convert it back to a list of raster layers for Cell Statistics to actually process. As a result I get this ERROR 000732:

Script failed because:      ERROR 000732: Input Raster: Dataset Development;Vegetation;Hydrology;Elevation does not exist or is not supported at this location:

I've tried reading in 'input' with .values instead of .valueAsText, so it comes in as a list, but then I get the error that 'cannot concatenate 'str' and 'list' objects' in the getattr line.

Any ideas what I can do about this?

I know this overall structure may seem unnecessarily complicated; but it's a key part of the bigger project. I'm not attached to specifically using getattr here, but solutions which preserve the basic structure, namely calling functions from a custom module using a string of the function's name, would be ideal


texts:
I need to convert DWG/shp files to some vector data file that is supported in
Openlayers (e.g KML,GML). I tried following this guide but it doesn't work.

Please advise me some converter (not necessarily free) that I can use.


texts:
I would like to explore outbreaks of Legionnaires disease in relation to cooling tower locations. This work would be confined to Colorado, potentially even just the Denver metro area.

I am wondering if anyone has ideas on where I might find data on cooling tower locations for Colorado?

If no dataset currently exists, I am considering digitizing cooling tower locations from satellite imagery, but that would pose obvious challenges.


texts: Downloading Sentinel 2 imagery using Semi-Automatic Classification Plugin?

texts: I have been following this video (https://www.youtube.com/watch?v=FcETq8OWM0k) to download Sentinel 2 imagery using the Semi-Automatic Classification Plugin in QGIS 3. However, when I click run and select the folder to download the files into, there is just a beeping sound and the files do not download. These are the settings I have been using:


texts: ArcFM AUs not firing when scripted from Python

texts:
I have a Python script to edit data that is "minerfied" using ArcFM.  The items to be edited require AUs to be fired off.
Example:   a field user uses ArcGIS Collector to edit some attributes on a point.  The Python script then takes the field edits, uses the special code from Schneider to enable ArcFM editing in my script, then make the edits to the minerfied data.
In the sample script Schneider has provided, there is a line that says:
au.AutoUpdaterMode = mmAutoUpdaterMode.mmAUMNoEvents

I believe the mmAUMNoEvents part tells the script to NOT fire AUs.  A colleague suggests that the last part of that line should instead be mmAUMStandAlone.  When I do this, my script crashes when it attempts to edit the data (the actual arcpy function fails).  I get this:
File "D:\Workspace\XXXXX\XXXXX.py", line 198, in <module>
uCursor.updateRow(urow)

File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\arcobjects\arcobjects.py", line 102, in updateRow
    return convertArcObjectToPythonObject(self._arc_object.UpdateRow(*gp_fixargs(args)))

Does anyone know of a good resource online that discusses the mmAutoUpdaterMode item in detail?  So far, I can't find much.
If anyone else has had this issue, or a similar issue, please let me know what you did to troubleshoot this and if you came up with a solution.


texts: Inserting full bleed images in browser background instead of maps using CartoDB?

texts:
I am building a product that will utilize high fidelity photography.

I would like to show my photos in a slideshow format that takes up the entire area of the browser.

CartoDB and the Odysssey.js would work well.

Can I substitute the maps for photography in the background of my browser?


texts: Joining intersected polylines using rules based rendering

texts:
Possible Duplicate:
  Setting rendering order of line layers in QGIS

I have trouble getting two intersecting polylines to look like they were joined at rendering. I am using rules based rendering and these two lines are rendered with same rules, so they look identical. One line ends where the other begings, but the line does not render to look smooth - rather it renders the end of other line.

Example:


texts: QgsExpression.evaluate() is deprecated from PyQGIS?

texts:
I am developing a plugin for QGIS 2.18 and I am sending queries to a layer:

self.userLayer.startEditing()

for feature in self.userLayer.getFeatures():
    value = expression.evaluate(feature)
    feature["ID"] = value
    self.userLayer.updateFeature(feature)

self.userLayer.commitChanges()

and I get a DeprecationWarning:

DeprecationWarning: QgsExpression.evaluate() is deprecated

What is the currently preferred way to get an expression evaluated?


texts: Looking for vb.net code to read and show attributetable fileld alias names in arcgis

texts:
have anybody codesnippets or complete code in vb.net to read and show the alias fieldnames of an attribute table of an shapefile in arcgis?

To open the Shapefile, I have the following code:

Public Function PersonalGeodatabaseÖffnen() As IWorkspace
        Const path As String = "C:\\DatabaseGIS\xxx.mdb"
        Dim workspaceFactory As IWorkspaceFactory2 = New       SapefileWorkspaceFactoryClass
        Dim workspace As IWorkspace = workspaceFactory.OpenFromFile(path, 0)
        Return workspace
        MsgBox(path)
    End Function


texts: ArcMap error drawing postgis layer - The number of points is less than required for feature

texts:
I have one big polygon called huc and a multilinestring called streams, I am running an intersection on between the two layers to clip the streams to the huc.

here is what I have so far

drop table if exists base_layers.streams_huc_clip;
create table base_layers.streams_huc_clip as
select s.*,case
    when st_within(s.shape,huc.shape) then s.shape
    when st_intersects(s.shape,huc.shape) then st_intersection(s.shape,huc.shape)
    end as geom2
        from postgres.streams s join base_layers.huc_hl_dissolve huc on st_intersects(s.shape,huc.shape);

--pretty straightforward query in postgis, this query returns with no errors and correctly gets displayed in QGIS.

However when I try and bring it into ArcGIS I get this error and it does not draw

One or more layers failed to draw:

lucz_2017.base_layers.streams_huc_clip:  The number of points is less than required for feature

so here are the following steps I tried after to do after it would not show up in ArcGIS to no avail.

dropped the original geometry so I just had one geometry column, renamed the new geom to shape

alter table base_layers.streams_huc_clip drop column shape;
alter table base_layers.streams_huc_clip rename geom2 to shape;

dropped the objectid and added new primary key, changed owners,

alter table base_layers.streams_huc_clip drop column objectid;
ALTER TABLE base_layers.streams_huc_clip ADD COLUMN id SERIAL PRIMARY KEY;
alter table base_layers.streams_huc_clip owner to base_layers;

created an index and tried updating the srid. none of this worked!

ALTER TABLE base_layers.streams_huc_clip ALTER COLUMN shape type geometry(multilinestring, 3424) using st_multi(shape);
CREATE INDEX streams ON  base_layers.streams_huc_clip USING GIST (shape);

also ran this query

SELECT * FROM public.geometry_columns ORDER BY type

which returned

I am thinking about running an st_dump on the geometry to turn it into a linestring, I will take any suggestions on why this layer is not displaying and what I can do to make it work

I ran an st_dump and turned the geometry into a linestring but I still got the same arcmap error

drop table if exists base_layers.streams_huc_clip;
create table base_layers.streams_huc_clip as
select (st_dump(geom)).geom as geom2,t.* from(
select s.*,case
    when st_within(s.shape,huc.shape) then s.shape
    when st_intersects(s.shape,huc.shape) then st_intersection(s.shape,huc.shape)
    end as geom
        from postgres.streams s join base_layers.huc_hl_dissolve huc on st_intersects(s.shape,huc.shape))t;

ST_MakeValid() did not help either

I make edits to these tables outside the ESRI stack.


texts: Is it possible to export an image that can be zoomed-in on?

texts: I am using QGIS on a Mac.  I would like to export my map so that the viewer can zoom in on it, just as I view the map on my QGIS canvas.  I am not sure how to word what I want.  If you see the map from the zoomed-out perspective, the symbols all bunch together, you need to be able to zoom in to see the details.  Anyone have any idea if this is possible?  I tried to export a high res image in the proint composer, but its not nearly enough resolution.


texts: Algorithm to find out points of inflection for a polyline

texts:
I am trying to find out the points of inflection , i.e. points where the curves in a line start and end. If you look at the image, the green line may be a road or a stream, and the black points are the points where the curves start and end.

What would be the high level steps to automate the generation of these points? I have ArcGIS desktop and am quite handy with ArcObjects.


texts:
I am trying to calculate the shortest path's distance from certain locations to their closest schools,
I used the Pgr_WithPointsCost Function to get the result from one location to its 10 closest destinations (One to Many), But I can't use the (many to many) scripts since Pgr_WithPointsCost will calculate the shortest distance between every point of the source and every target point. The problem is my data set is huge and when I launch the script it crashes.
This is the script I used to calculate the distance


texts: Determining rotation of rectangle features for Data Driven Pages in ArcGIS Desktop?

texts:
I am trying to generate a large number of plans using the Data Driven Pages in ArcGIS Desktop.

I have created an index layer comprising a series of rectangle features that correspond to set scales on my layout page (eg 1:100, 1:500, etc).  Due to the odd arrangement of the features I laid out the plans with rectangles that are often rotated.  There wasn't a particular feature that I could use alternatively to create a nice layout using strip or grid tools.

I can see in the data driven page there is an option to include rotation.  However I cannot figure out how to measure the rotation of the rectangle features.

Is there a way to do this automatically in the attribute table, similar to the way you can calculate geometry?


texts:
I have a question about the code below. Although I've defined geometry as a triangle, when I download the file, the image has a square shape. I running that code using node js. Can you help me, please?

  ee.data.authenticateViaPrivateKey(config.gee, () => {

    ee.initialize();

    var geometry = ee.Geometry.Polygon([[[-58.22065005792236, -38.39475941140011],[-58.21198115838623, -38.391698569574615],[-58.222109179626464, -38.38813303001853]]])

    var collection = ee.ImageCollection('COPERNICUS/S2').filterBounds(geometry).filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 5));

    var cloudMasked = collection.map(image => {
      var cloud = ee.Number(2).pow(10).int();
      var cirrus = ee.Number(2).pow(11).int();
      var qa = image.select('QA60');
      var mask = qa.bitwiseAnd(cloud).eq(0).and(qa.bitwiseAnd(cirrus).eq(0));
      return image.updateMask(mask);
    });

    var ndvi = cloudMasked.map(image => {
      return image.addBands(image.normalizedDifference(['B8', 'B4']).rename('NDVI'));
    });

    var reducer = ee.Reducer.percentile([90]);
    var ndvi90 = ndvi.select(['NDVI']).reduce(reducer).clip(geometry);
    var url = ndvi90.getDownloadURL({'region':geometry,'scale': 10})

    console.log(url);
  });


texts:
I'm using WFS protocol on Vector layer / GeoServer / PostGIS

Currently the map is refreshing and calling new requests whenever i do anything like searching or drawing and when i move back to a part of the map which previously had features already downloaded, it refreshes and call new request again. i need to  reuse previous data for previous requests.

I've read that the OpenLayers does not have this functionality. but this was a long time ago.

Is there a way to cache wfs returns by check the boolean in top of the code or anything like that ?

That will help me implement my BACK button so that i can store previous requests in a stack and pull requests when the user clicks the back button.


texts: Remediating path that crosses international dateline

texts:
I'm trying to plot the paths of multiple ships on a world map using R and a series of Lat/Lon points.

All is well until the ship crosses the international dateline (-180/180) when the path jumps across the map.

I tried applying the st_wrap_dateline() function in R but it looks like it only works on datasets with two points, a start & end coordinate.

My R code is below as well as a screenshot of the plot produced.

# Download background Blue Marble Globe image
download.file("https://www.researchvessels.org/images/nasa_base_v2.png", "nasa_base_v2.png", mode="wb")
# Download ship lat/lon data
download.file("https://www.researchvessels.org/images/shipdata.RDATA", "shipdata.RDATA", mode="wb")
load("shipdata.RDATA", envir = parent.frame(), verbose=TRUE)

# Ensure that the required packages are installed
list.of.packages <- c("ggplot2", "sf", "dplyr", "png", "grid")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)>0) {install.packages(new.packages)}

library(ggplot2)
library(sf)
library(dplyr)
library(png)
library(grid)

bluemarble_bg <- png::readPNG("nasa_base_v2.png")

xlim = c(-180.0,180.0)
ylim = c(-90.0,90.0)

shiptracks$Date.Time <- as.POSIXct(shiptracks$Date.Time, tz="GMT", origin="1970-01-01")

pal <- c("ShipA" = "#488f31",
         "ShipB" = "#FF00FF",
         "ShipC" = "#fff1a9",
         "ShipD" = "#f19d61",
         "ShipE" = "#de425b")

shiptracks %>% group_by(Vessels.Name)

shiptracks <- st_as_sf(shiptracks, coords=c("Lon", "Lat")) %>% st_set_crs(4326)

shiptracks <- cbind(shiptracks, st_coordinates(shiptracks))

shiptracks <- shiptracks %>%
  st_sf() %>%
  st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"),
                   quiet = FALSE) %>%
  sf::st_sf(crs = 4326)

shipmap <- ggplot(shiptracks, aes(x=X, y=Y, group = Vessels.Name)) +
  coord_sf(xlim = xlim, ylim = ylim, expand = FALSE,
           crs = 4326,
           datum = sf::st_crs(4326), label_graticule = waiver(),
           label_axes = waiver(), ndiscr = 100, default = FALSE,
           clip = "on") +
  annotation_custom(rasterGrob(bluemarble_bg,
                               width = unit(1,"npc"),
                               height = unit(1,"npc")),
                    -Inf, Inf, -Inf, Inf) +
  geom_path(data = shiptracks, group = shiptracks$Vessels.Name, color=pal[shiptracks$Vessels.Name],
            aes(x=X, y=Y), cex=1, show.legend = TRUE) +
  geom_point(data = shiptracks, group = shiptracks$Vessels.Name, color=pal[shiptracks$Vessels.Name], aes(x=X, y=Y), cex=1, show.legend = TRUE) +
    theme(axis.line = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.background = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "mm"),
        panel.spacing = unit(0, "mm"),
        legend.position = c(-120, -20),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 5),
        legend.key.size = unit(0.01, "npc"))

shipmap


texts: How to find and delete all road unattached road segments

texts:
I have a large road network that I'm using for routing and agent-based modelling and I've just noticed that it has a lot (possibly hundreds) of road segments that are not connected to the full network. I think this was because the road network was clipped to fit a county line that was a bit wiggly. My agent-based model places agents on road segments at random and when it places them on unattached road segments it fails.

Is there a way to highlight and remove all unattached road segments using QGIS or ArcGIS?

I've done some Googling and found info on 'dangles' but I've not been able to get the job done.


texts: How to split the polygon, operation using of JSTS?

texts:
I'm use JSTS 0.17.0 (https://github.com/bjornharrtell/jsts).
How to split the polygon (polygon to multipolygon or two polygons) ?

Input:

POLYGON ((1 1, 1 9, 9 9, 9 1, 1 1))
LINESTRING (0.5 5, 9.6 5, 9.6 3.9)

Result (or):

MULTIPOLYGON (((9 5, 9 1, 1 1, 1 5, 9 5)),((9 9, 9 5, 1 5, 1 9, 9 9)))
GEOMETRYCOLLECTION ( POLYGON ((9 5, 9 1, 1 1, 1 5, 9 5)), POLYGON ((9 9, 9 5, 1 5, 1 9, 9 9)) )

JTS does so:

How do too JSTS ?

PS2
Sorry for my English.


texts: QGIS 3x Label doesn't move when moving the point feature

texts:
I want to move my label and basically set some expressions for it.
Unfortunately since I started to use the data attribute table auxiliary columns for location and movement, my layer is not really editable this way.
First of all I can't move the selected point layer, to which the label is assigned, like below (see the yellow dot).

Secondly the label placement (offset) doesn't work anymore.
In turn, my street name shifts from one side to another as the scale changes.

All expressions regarding to font size algorithms, rotation etc stopped working beyond this point.
Is there any way to combine the auxiliary rotation with rotation and offset set by expression?


texts: Consistent layer colors that are automatically adjusted based on adjacent polygons?

texts:
I have several polygon layers and I want to set 4 colors to fill all of those, but the color of one layer should not be the same with the adjacent polygon.

I have read several questions about randomizing colors or using color ramps but most appear to be reliant on attribute fields/columns. It would be very tedious to make another style column in the attribute table, and in my case, it would be impractical because I want the colors to automatically adjust so that they are distinct from each neighboring polygon.

Is there a way to do this?


texts: How to convert geometry of shapefile to decimals

texts:
I try to create a line with two points. These points are obtained in a shapefile layer.

How can I convert that points to decimals?


texts: Porting ArcGIS Engine to ArcGISRuntime SDK for .NET, dealing with interior/exterior polygon rings

texts:
I've inherited some code that converts ArcGIS Engine IPolygon4 objects to GeoPackage-compatible objects. The format is very specific in terms of what's permitted with exterior and interior rings.

I need to port this code to ArcGIS Runtime for use on Android. The API for the Runtime polygon class Polygon is really sparse, and contains no properties or methods related to the polygon's rings. The only property I have is an IReadOnlyList<ReadOnlyPart> holding the polygon's parts.

I don't know if this collection is ordered, or if I can order it in some way to determine which ring is the exterior one.

Is there any way that I can identify the interior/exterior rings of a Polygon in ArcGISRuntime SDK for .NET?

Here is the code I am working with. WkbMultiPolygon is the custom polygon object I am trying to convert to.

using ESRI.ArcGIS.Geometry;
using System.Linq;

public static WkbMultiPolygon CreateMultiPolygon(IPolygon4 esriPolygon) {

    var exteriorRingGeometryCollection = (IGeometryCollection)esriPolygon.ExteriorRingBag;

    var exteriorRingCount = exteriorRingGeometryCollection.GeometryCount;

    // create my list of custom polygons
    var wkbPolygons = new List<WkbPolygon>(exteriorRingCount);

    for(var exteriorRingIndex = 0; exteriorRingIndex < exteriorRingCount; ++exteriorRingIndex) {
        var exteriorRingGeometry = exteriorRingGeometryCollection.Geometry[exteriorRingIndex];

        var interiorRingCollection = (IGeometryCollection)esriPolygon.InteriorRingBag[exteriorRingGeometry as IRing];

        // create my custom polygon object and add to list.
        wkbPolygons.Add(new WkbPolygon(ToLinearRing(exteriorRingGeometry as IPointCollection),
                                       Enumerable.Range(0, interiorRingCollection.GeometryCount)
                                                 .Select(index => interiorRingCollection.Geometry[index] as IPointCollection)
                                                 .Select(ToLinearRing)
                                                 .ToList()));
    }

    return new WkbMultiPolygon(wkbPolygons);
}


texts:
I am working with 2 shapefiles. One contains the boundaries of a city, the other contains polylines representing the topographic features of the area. Essentially, the polylines make a topographic map.

I want to calculate the average elevation of the city. I have some basic GIS knowledge, so I can select out just the polylines in the city... but I'm not sure how to determine the city's average elevation using this data? I could average the polylines together controlling for their length, but that's not quite accurate - a long line that goes around a lake contains a lot more area at the same elevation than a long line that goes around a big hill, for example.

Suggestions on how to do this?


texts: Why are the labels not displaying the correct (virtual) attribute values?

texts:
I have a polygon and I create a new field area_ha using this expression in the field calculator:
round( $area /10000,1)
When I use that field for my labels, everything is fine and works as it should work:

I create a new field, but this time a virtual field with the name "area_ha_virtual_field". I use the same expression in the field calculator. The resulting value in the attribute table is fine. But when I want to use that new virtual field to label my blue polygon, it shows a totally different number than in my attribute table:

Also, I get that same weird number as a label, when I use that expression in my labeling properties:

My CRS is 25832.


texts: Calculating distance between point (address) and every road in city using QGIS

texts:
I'm new on QGIS, and I'm trying to measure the distance between an address and every road in the city.

I used NNjoin but I got only the shortest distance between the address and the nearest road.

Is there a package or plug-in that I can use to calculate a distance between the address and all roads in the city?


texts: Using QGIS to remove portions of polygons that overlap in same layer?

texts:
Is there a way to easily eliminate the overlapping portions of polygons in the same layer?

For example, I want to turn this

into this

I used Union to break the polygons apart where they overlapped, and then used Delete Duplicate Geometries to get rid of half the overlaps, but still had to manually delete the others. That's fine for the example, but won't work with the large amount of polygons that will be in the real data.


texts:
I am looking to perform a simple (one observer location) viewshed analysis on a standard DEM using GDAL (either python or .NET bindings).

I have searched, but haven't found any readily available implementation.  Can anyone point me in the right direction?


texts:
I tried to get WMS tiles from geoserver like

/geoserver/gwc/service/tms/1.0.0/layername@EPSG%3A4326@png/15/19142/9530.png

or with another projection

/geoserver/gwc/service/tms/1.0.0/layername@EPSG%3A900913@png/15/19142/9530.png

But all the time I got empty tile. For grabbing them I use openlayers.
All in all the tile, which should be in 19142/9530.png are located in:

geoserver/gwc/service/tms/1.0.0/layername@EPSG%3A900913@png/15/19801/22523.png

So as I could understood all the tiles have pretty big offset

The offset tiles location

and I got them in the center of the map.That layer which should be in Europe

Where tile are really located

are located in the part of Africa and ocean(link 1).

That's why I  couldn't use them with OSM maps and ect(because my tiles are simply empty).

Can someone help me to manage with this problem? Is this a issue of geowebcache (geoserver) or I'm doing something wrong?


texts:
I am not sure what happened between where I was working last week and what I am seeing this week. It may be that QGIS updated without my knowledge, but I am not sure because I think the version is still the same though I have seen some things that I thought look different.

I had enabled the OpenLayers Plugin (using Google Streets) and was using a scale of 1:186. There is also a georeferenced image file layer over the Google street layer. I am now unable to see the Google street layer at the scale I had been mapping streets and other points from my image file (there is enough overlap to track changes to the city I am studying). In addition, everything seems skewed slightly to the left, or up, or down from where I had it at the 1:186 scale.

I am not sure what the issue is I need to fix here. Does anyone have any advice or tips?


texts: How to install OpenLayers plugin in QGIS1.8

texts: I am using QGIS 1.8.0 on Windows 7 (32-bit) but am unable to install the OpenLayers plugin.  The Python plugins menu does not include OpenLayers on the plugin list.  Can it be somehow added manually?


texts: ogr2ogr: Create CSV with strings wrapped in double quotes

texts: Is there a way to force ogr2ogr to wrap string-column values with double-quotes?  I have a column with string values like "000010001", which OGR correctly interprets as a string when writing to CSV.  However, Excel converts that value to a numeric 10001 when parsing the CSV, due to a lack double-quotes wrapping the value.


texts: Could use a bit of XAML guidance when designing with ESRI Silverlight Viewer

texts:
I'm trying to implement a version of this in my web map:

As you can see, part of the XAML is:

    <esri:Map x:Name="MyMap" WrapAround="True" ExtentChanged="MyMap_ExtentChanged" ExtentChanging="MyMap_ExtentChanged">
        <esri:ArcGISTiledMapServiceLayer ID="StreetMapLayer"
            Url="http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"/>
    </esri:Map>

For me, the problem lies in that I'm trying to create a new widget that will show this information in a little window that pops up when someone opens it. The layers are defined elsewhere in the project (in an AppConfig.xml file), I don't need to create a new map in this widget. I just want something that can get the extents that the user is zoomed into. That's what I don't know how to do... I'm not sure how to get that information if I'm not defining the map layers in that particular XAML file I created for the widget.

I apologize if there's some really easy answer. I've been working on this web map stuff for only a couple weeks and I run into this type of issue a lot, where sample code is shown defining layers and going from there... but it's not necessarily transferable to designing with widgets. Or maybe it is, and I'm just a complete newbie to Silverlight, C# and XAML :(.

Oh, and I'm working on adding a widget to this.  It'll just go next to the Print one =).


texts:
I have two layers in my Project

1)Points Layers

2) Polygon Layer

I need to apply nearest neighbour Analysis between point and Polygon, but i am unable to get .

I tried the below code.I am using Qgis 3.0

from mmqgis import mmqgis_library as mmqgis

As soon i added this line in my python console it got

  No module name mmqgis found


texts:
For QGIS 2.10.1 Pisa, Win8

I am attempting to plot my data (about 50 points, imported CSV, converted to shapefile); as graduated by size, divided by natural breaks, 8 classes, link class boundaries, simple circle markers.

The problem is one point will not plot, unless the class is only that point (class 1):

1     .449872-.449872
2     .449872-.649973
3     .648873-.749909
etc.

The point should plot with class 2 but does not. If I delete class 1, the point goes away.  I have tried tweaking the class values such that the value should be included in class 2 also, but no difference.

The circles are sized with the lowest values as the largest circles, if that matters. Additionally, sometimes the “show feature count” function gives the correct number of feature in the class but other times, it only lists zeros, (not a separate question, maybe a related symptom).


texts: I have problem editing data in ArcPad because of the limitation of software. Can you help me? My client wants to do editing in ArcPad but because of the limitation, I am trying to find a solution to help them by using AXF file out of ArcPad. Do you know how to edit axf file in ArcMap ? But I am not working in shapefile but checkout data from feature layer which is carrying subtype and domain. Looking some help from you.


texts: Using raster calculator to create cumulative ZTV in QGIS

texts:
I need to create a raster layer with cumulative values for 72 viewsheds.
Each viewshed is from a different point along a line and has been clipped to a 20km buffer from each viewing point.
When I try to add them all together in raster calculator in QGIS, I am left with a raster that only covers where each viewshed overlaps (see image below with 3 of the rasters covering different areas and the result from the raster calculator in red where they overlap.).
I need the coverage to be over the whole area so it will tell me where only 1 point is visible in an area all the way up to potentially all 72 being visible.
Is there a way to do this in QGIS?


texts:
I want to load an OSM map via de XYZ Tiles layer, but is is not showing.
And I can not figure out what it is. What can I check to see what is wrong?


texts: Enabling snapping for OnLine method of ArcPy Python Add-In?

texts:
I am working on a Python Add-In using ArcGIS 10.2.2 where the user can draw lines and the Add-In will add dimension annotation showing length for each segment of the line that is drawn.  I can get this working just fine and the text is angled properly, but the OnLine() method does not use snapping in the map.  The add-in does however, properly set the snapping (I verified the snapping tolerance in the map and tested with different tolerances) in the map for everything BUT the onLine() method for the tool.

In the ArcObjects help docs, it says this is usually invoked on a OnClick method from the ICommand Interface.  I do not know how I could use this in a Python Add-In.  Here is my attempt to use comtypes and some custom modules to invoke snapping for my Add-In.  Unfortunately, it is not working.  Maybe I should attempt this in .NET?

import arcobjects #custom module
from arcobjects import elements #custom module
import arcpy
import pythonaddins
import math
import comtypes.gen.esriControls as esriControls
import comtypes.gen.esriArcMapUI as esriArcMapUI
from arcobjects import CType, NewObj, CLSID #custom

def midpoint(point_a, point_b):
    x1 = point_a.X
    y1 = point_a.Y
    x2 = point_b.X
    y2 = point_b.Y

    # Find midpoint
    x = (x1 + x2) / 2.0
    y = (y1 + y2) / 2.0
    return (x, y)

def get_angle(xy1, xy2):
    """Calculate azimuth angle from two points. (Zero is north.)"""
    # ArcPy point objects
    x1, y1, x2, y2 = xy1.X, xy1.Y, xy2.X, xy2.Y
    dx, dy = (x2 - x1, y2 - y1)
    return math.degrees(math.atan2(dy, dx))

class DrawLine(object):
    """Implementation for AddMeasurementAnno_addin.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.shape = 'Line'
        self.cursor = 3

        # arcobjects get current app (IApplication interface Pointer)
        self.pApp = arcobjects.GetCurrentApp()
        self.pDoc = self.pApp.Document
        self.pMxDoc = CType(self.pDoc, esriArcMapUI.IMxDocument)
        m_SnappingEnv = NewObj(esriControls.Snapping,
                               esriControls.ISnappingEnvironment)
        m_SnappingEnv.Enabled = True
        m_SnappingEnv.IgnoreIMSLayers = True
        m_SnappingEnv.ShowSnapTips = True
        m_SnappingEnv.ShowSnapTipBackground = True
        m_SnappingEnv.Tolerance = 10
        m_SnappingEnv.SnappingType = esriControls.esriSnappingTypeVertex + \
                                     esriControls.esriSnappingTypeEndpoint + \
                                     esriControls.esriSnappingTypeEdge
        print m_SnappingEnv.SnappingType
        m_Snapper = m_SnappingEnv.PointSnapper
        m_SnappingFeedback = NewObj(esriControls.SnappingFeedback,
                                    esriControls.ISnappingFeedback)
        m_SnappingFeedback.Initialize(self.pMxDoc, m_SnappingEnv, True)

    def onLine(self, line_geometry):
        # Report back mouse position at line endpoint
        #
        m_CurrentMouseCoords = self.pMxDoc.CurrentLocation
        print m_CurrentMouseCoords.X, m_CurrentMouseCoords.Y

        # loop through geom and grab vertices
        for part in line_geometry:
            for i in range(part.count):
                if i != (part.count - 1):
                    pt1, pt2 = part.getObject(i), part.getObject(i+1)
                    line = arcpy.Polyline(arcpy.Array([pt1, pt2]))
                    length = str(round(line.length,1)) + "'"
                    x, y = midpoint(pt1, pt2)
                    elements.add_text(self.pApp, text=length, name='length', size=9, bold=True,
                                      x=x, y=y, angle=get_angle(pt1, pt2), view='data')

Here is what it does in ArcMap, but the distances are not correct because the line was not snapped.  Without snapping this Add-In is pretty useless:

And the snap options window:

Has anyone had success in doing this?


texts: I have a .sid image but can't figure out the CRS

texts:
There is a .map file:

OziExplorer Map Data File Version 2.2
Tiritir_2.sid
C:\OziExplorer\Maps\Islands\Tiritir_2.sid
1 ,Map Code,
Geodetic Datum 1949,WGS 84,   0.0000,   0.0000,WGS 84
Reserved 1
Reserved 2
Magnetic Variation,,,E
Map Projection,(NZG) New Zealand Grid,PolyCal,No,AutoCalOnly,No,BSBUseWPX,No

my best guess is NZGD49 North Island but when I use this projection the data from the GPS does not match.  Any ideas?


texts: QGIS, clip my vector data with a DEM raster

texts:
I know that I cannot simply clip vector data with raster data.

I'm working in QGIS and I have vector point data of earthquakes (XYZ), and I have it projected on top of my raster data (.tif). However, my point data goes way beyond my raster data boundaries which is unnecessary.

How can I clip out the vector point data that is beyond the boundaries of my DEM file?
They are 100% in the same CRS.  I have tried "polygonize" (raster to vector) but for some reason, it kept giving me errors saying that the file couldn't be created. I also tried vector>research tools>select by location, which didn't work either.


texts: Grass algorithm r.neighbors doesn't work in Qgis 3.4. How can I solve this?

texts:
When I try to run the algorithm it gives me this error:

  Traceback (most recent call last): File
  "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py",
  line 411, in processAlgorithm Grass7Utils.executeGrass(self.commands,
  feedback, self.outputCommands) File
  "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py",
  line 372, in executeGrass for line in iter(proc.stdout.readline, ''):
  File "C:\OSGEO4~1\apps\Python37\lib\codecs.py", line 322, in decode
  (result, consumed) = self._buffer_decode(data, self.errors, final)

  UnicodeDecodeError: 'utf-8' codec can't decode byte 0xe8 in position 50: invalid continuation byte


texts: Merging all feature classes with same name nested in multiple geodatabases using ArcPy says feature class already exists?

texts:
I'm following very closely this topic Merging all feature classes with same name nested in multiple geodatabases using ArcPy?, it seems very useful to my python script

this is the script I following from made by Emil Brundage (answered Sep 8 '17 at 15:21):

import arcpy
import os

workspace = "C:/Wiley/P1/gis"
outdir = "C:/Wiley/P1/gis/HK80.gdb"

walk = arcpy.da.Walk(workspace, datatype="FeatureClass", type="Polygon")
di = {}

for dirpath, dirnames, filenames in walk:
          for filename in filenames:
             fcFullPath = os.path.join(dirpath, filename)
             #try adding feature class to dictionary with already-existing key
             try: di [filename] += [fcFullPath]
             #key not in dictionary yet. Create key
             except: di [filename] = [fcFullPath]

#iterate dictionary
for filename in di:
    #get merge feature classes
    mergeFcs = di [filename]
    output = os.path.join(outdir, os.path.basename(fcFullPath) + "_merge")
    arcpy.Merge_management(mergeFcs, output)

However, when I run the script I get a message error show up, a feature class already exists. The script does exactly what I want, but how could I fix the script to avoid this error type.

Following the suggestion from amasephy add the arcpy.Exists and arcpy.Delete_management function on my script:

workspace = r'...\GDBESTATALES_PRUEBA'
outdir = r'...\GDB_Nacional.gdb'

walk = arcpy.da.Walk(workspace, datatype="FeatureClass", type="Polygon")
di = {}

for dirpath, dirnames, filenames in walk:
          for filename in filenames:
             fcFullPath = os.path.join(dirpath, filename)
             #try adding feature class to dictionary with already-existing key
             try: di [filename] += [fcFullPath]
             #key not in dictionary yet. Create key
             except: di [filename] = [fcFullPath]

#iterate dictionary
for filename in di:
    #get merge feature classes
    mergeFcs = di [filename]
    output = os.path.join(outdir, os.path.basename(fcFullPath) + "_merge")
    arcpy.Merge_management(mergeFcs, output)
    if arcpy.Exists(output):
     arcpy.Delete_management(output)

However, nothing happens what I want it, the only one happens is the script takes the last geodatabase and save the last feature class in the other feature class (GDB_Nacional.gdb), that's all.


texts: Overwriting files in geodatabase outside workspace environment using ArcPy?

texts:
I initially set the workspace environment to a directory and included the option, arcpy.env.overwriteOutput=True, but at some point in my script, I convert a shapefile in that directory to a geodatabase file. When I do this, I get the error:

  File (filename) already exists.

I assume that because the overwriteOutput option is specified for the workspace, it doesn't apply to the geodatabase within that directory. However, I also tried to change the environment to the geodatabase and set the overwrite option to true again. The error remained.

How can I fix this (without deleting the file)?


texts: Merge lines while adding attributes spatially

texts:
What I have:

A shapefile with multiple lines.

Example:

line 1: A -> B -> C with attributeX=5

line 2: A -> B -> D with attributeX=7

What I want

Lines between each nodes with added attributes.

Example:

A -> B with attributeX=12

B -> C with attributeX=5

B -> D with attributeX=7

(QGIS 2.0)


texts: Loading an AGOL map with a token in ArcGIS 100.*

texts:
I am trying to upgrade from ArcGIS 10.* to the new 100.* framework and I am having trouble getting the secured map to load. This is the old code:
_credential = [[AGSCredential alloc]initWithToken:token];
if (_credential != nil)
{
    self.webMap = [[AGSWebMap alloc]initWithItemId:mapId credential:_credential];
    self.webMap.delegate = self;
    [self.webMap openIntoMapView:self.mapView];
}

Here is what I am trying currently:
_credential = [[AGSCredential alloc]initWithToken:token referer:nil];
if (_credential != nil)
{
    AGSPortal* portal = [AGSPortal ArcGISOnlineWithLoginRequired:false];
    portal.credential = _credential;
    [portal loadWithCompletion:^(NSError * _Nullable error) {
        if(error)
        {
            NSLog(error.localizedDescription);
        }
    }];
    AGSPortalItem* portalItem = [[AGSPortalItem alloc]initWithPortal:portal itemID:mapId];
    self.webMap = [[AGSMap alloc]initWithItem:portalItem];
    self.mapView.map = self.webMap;
}

The problem is if I set the login required to false I get a 400 not logged in error. If I set it to true it pops up a login alert. I think I'm not creating the credential correctly maybe? Like I don't know what to put for the referer specifically. The docs on this are sparse.


texts:
I have a centre point XY. I have a series of measurements around the point in xy offsets in centemetres from XY. I need to create a polygon of these measurements.

With the first measurement x is always = 0. Sometimes the first measurement is not due north from XY. When this is the case I have an angle describing north in relation to the first measurement, e.g. 325, if the first measurement is on an azimuth 35 degrees east of north.

How do I adjust x and y to get their "true" positions in relation to XY?

I've tried (in python):

angle =  math.radians( 360 - 325 )
newX = X + x * cos(angle)
newY = Y + y * sin(angle)

but things don't quite look right.


texts: Adding image links in GeoJSON object to get images in popup using OpenLayers 3?

texts:
I have added a geojson file on my map. By clicking on a geojson feature I am getting a pop up which is displaying geometry and coordinate values. Now I want to display image in a pop up also.

How to add it?

This is the geojson code:

{
  'type': 'FeatureCollection',
    'crs': {
      'type': 'name',
        'properties': {
          'name': 'EPSG:4326'
        }
    },
      'features': [{
        'type': 'Featurejson',

        'geometry': {
          'type': 'Point',
          'coordinates': [92.460416718365195, 10.845240354970281],
          //'imgurl': 'http://www.w3schools.com/html/default.asp'
          //'name':"file:///C:/Documents%20and%20Settings/sonal/Desktop/html/point.PNG"
        },
        "properties": {
          //"image": "http://www.w3schools.com/html/default.asp",
          //"href": "http://www.w3schools.com/html/html5.gif",
          "name": "file:///C:/Documents%20and%20Settings/sonal/Desktop/html/point.PNG",
          //"1":"point.png",
          "description": "43.732, -83.9387"

        }
      },


texts: I'm managing a project - http://nepanode.anl.gov aimed at improving data awareness and sharing around Environmental impact analysis. I'm struggling to come up with ways to engage users and show them how they can apply the site. I'm looking for specific ideas, games, demos, tutorials to demonstrate how the site can be used; looking for new ideas and/or links to existing training resources that aim to engage users through games or case studies? Limited to users only using the site and/or using QGIS


texts: Lengthening 1st and last point of GeoJSON LineString containing lat/lon coordinates maintaining angle and direction

texts:
I am fairly new to the gis topic and may not be using correct terms.

I have a GeoJSON LineString containing lat/lon coordinates. I would like to lengthen the 1st coordinate and the last coordinate by a certain distance, maintaining the angle and direction (heading?) between the 1st and 2nd coordinate respectively the "next to last" and last element.

in other words, lengthen the line between the 1st and 2nd coordinate respectively lengthen the line between the "next to last" and last coordinate. as a result, the coordinate of the 1st and last element changes.

The following is a simplified LineString:

from shapely.geometry import LineString

line_orig=LineString([(4,5),(6,7),(8,7),(10,8),(12,8),(13,7)])

The result would be this one if the ends are lengthened by 1:

line_edited=LineString([(3.29,4.29),(6,7),(8,7),(10,8),(12,8),(13.71,7.29)])

not sure how to approach this. I would prefer a python library handling this server side. so far I was not successful finding something.

I would like to avoid calculating "manually" like

import math

# first line coords
firstLine_x1 = 4
firstLine_y1 = 5
firstLine_x2 = 6
firstLine_y2 = 7

# last line coords
lastLine_x1 = 12
lastLine_y1 = 8
lastLine_x2 = 13
lastLine_y2 = 7

def lineDirection(start_x1, start_y1, end_x2, end_y2):
    return math.atan2(end_y2 - start_y1, end_x2 - start_x1) * 180 / math.pi

firstLine_angle = lineDirection(firstLine_x2, firstLine_y2, firstLine_x1, firstLine_y1)
lastLine_angle = lineDirection(lastLine_x1, lastLine_y1, lastLine_x2, lastLine_y2)

print(firstLine_angle, lastLine_angle)


texts:
I have a C++ program running some OpenCV processing on 4096x4096 JPEG2000 images that have geographical metadata.  I don't want to modify the original, the results of my image manipulations are stored in a new image file.  Since my resulting image is missing the data, I have implemented a copy method in my program using CreateCopy() to copy the source image data into the new image.  On a 216kB image it takes 97sec to copy!

I have tried GDALSetCacheMax(200) and it seems to be worse.

I can't understand why such a small piece of metadata takes longer than even my OpenCV image processing algorithms!  Please help!

Here is a sample GDALINFO output for a correct image.  If I don't use CreateCopy, then my new image is missing the coordinate system, and the corner coordinates are just the pixel dimensions.

Driver: JP2OpenJPEG/JPEG-2000 driver based on OpenJPEG library
Files: <filename>.jp2
Size is 4096, 4096
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-97.025756835937500,32.398681640625000)
Pixel Size = (0.000001341104507,-0.000001341104507)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -97.0257568,  32.3986816) ( 97d 1'32.72"W, 32d23'55.25"N)
Lower Left  ( -97.0257568,  32.3931885) ( 97d 1'32.72"W, 32d23'35.48"N)
Upper Right ( -97.0202637,  32.3986816) ( 97d 1'12.95"W, 32d23'55.25"N)
Lower Right ( -97.0202637,  32.3931885) ( 97d 1'12.95"W, 32d23'35.48"N)
Center      ( -97.0230103,  32.3959351) ( 97d 1'22.84"W, 32d23'45.37"N)
Band 1 Block=1024x1024 Type=Byte, ColorInterp=Red
  Overviews: 2048x2048, 1024x1024, 512x512, 256x256
  Overviews: arbitrary
Band 2 Block=1024x1024 Type=Byte, ColorInterp=Green
  Overviews: 2048x2048, 1024x1024, 512x512, 256x256
  Overviews: arbitrary
Band 3 Block=1024x1024 Type=Byte, ColorInterp=Blue
  Overviews: 2048x2048, 1024x1024, 512x512, 256x256
  Overviews: arbitrary

GDAL code:

//GDAL Data
const char *pszFormat = "JPEG2000";
GDALDriver *poDriver;
GDALDataset *srcDataset;
GDALDataset *dstDataset;
GDALAllRegister();
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

//Get Source metadata
srcDataset = (GDALDataset *) GDALOpen(imagePath, GA_ReadOnly);

//Create new image data
imwrite(maskPath, dst);

//Write metadata to new image
dstDataset = poDriver->CreateCopy(maskPath, srcDataset, false, NULL, NULL, NULL);
if(dstDataset != NULL)
{
    GDALClose((GDALDatasetH)dstDataset);
}


texts: GME: could not open input point data source

texts:
I can not get this for loop to run in GME because
i is being changed to 1 and , out is not recognized:
wd <- "C:\Users\name\Desktop\KDE\2011\";
for (i in 1:15)
{
kde(in="C:\Users\name\Desktop\KDE\first15clip.shp"), out=paste(wd, "kde_p", 1, ".img"), bandwidth=1000, cellsize=10, where=paste("PName_ID=", 1));
}

Warning: An unexpected command argument was encountered: , out. The tool will attempt to run anyway, but may not produce the output you intended it to. Please check the results carefully.
WHERE clause: PName_ID=1
Error: The command text could not be interpreted. Please check the syntax of the command.
Error: An important error has occurred. Please include the information below if you submit a query about this error.
Length cannot be less than zero. Parameter name: length

Any ideas on how to fix this so the loop runs?


texts: GDAL SetGeoTransform does not work as expected

texts:
No transformation is visible after using SetGeoTransform to the dataset opened from a *.tiff file:

ds = gdal.Open(filename, GA_ReadOnly) # open the raster data set

geo_transform = (0.3,10,30,-10,80,-1) # complete random/arbitrary numbers
ds.SetGeoTransform(geo_transform)

cols = ds.RasterXSize
rows = ds.RasterYSize
band = ds.GetRasterBand(1)
data = band.ReadAsArray(0,0,cols, rows)

With plotting, I get the same results as plotting without the transformation:

plt.imshow ( data, interpolation='nearest', vmin=0, cmap=plt.cm.gist_earth)
plt.savefig('question.png',dpi=75)

GDAL version 1.11.3


texts: Subtract consecutive rasters from a list in Python

texts:
I have a list of 15 rasters as
b1.tiff
b2.tiff
b3.tiff
b4.tiff
I want to subtract the consecutive rasters, say b1.tiff from b2.tiff, b2.tiff from b3.tiff and so.
I see code with R but I want to do in Python. Can anyone help me with this?


texts: What is the time interval for daily GPM IMERG product? 0:00 to 23:59 local time? UTC?

texts:
I'm processing a bunch of GPM IMERG Daily products and comparing to historical station data. These old met stations report daily precip from 9am to 9am of the following day. So I'm wandering what is the time interval cor the daily GPM data? 0:00 to 23:59 local time?

Is there a service that will provide daily GPM data considering other time intervals (i.e. 9am to 9am following day)? Or will I need to download the 30min data and aggregate manually?


texts:
I am new to QGIS (ESRI shop at work).
I am attemping to concatenate fields in the label expression

field1 = number
field2 = number

label expression I tired but fails

str(field1) + " - " + str(field2)

it says its not valid. This would work in ArcMap... how is it different?


texts: How to apply trained classifier to entire Landsat collection in Google Earth Engine

texts:
I have classified the 2019 image using Google Earth Engine. How can I use the training data and code I have for 2019 and apply it from 1984 to 2018 instead of training/running a code for each image individually.

Link to my script: https://code.earthengine.google.com/15fc51687216911cfdf7d3c8c184938f


texts: How to access map service from GeoServer and can we use that service in ArcGIS API?

texts:
I have installed GeoServer and publish one shapefile on it. Now my question is, How can i access that shapefile from GeoServer? and can i work with that resource (shapefile) with ArcGIS API. As we can access any resource from ArcGIS sample server like the code below:

queryTask = new esri.arcgis.gmaps.QueryTask("sampleserver1.arcgisonline.com/ArcGIS/rest/services/…;);

Similarly, which piece of code i need to access my published file from GeoServer ?


texts:
I am new to coding and not too familiar with KMLs but I am trying to add an image icon overlay into an existing KML file (that contains polygon layer from shapefile) without having to overwrite the file. Previously Ive tried the "simplekml" python module which was giving me what I want but unfortunately I am not able to make  the modification to an existing KML file.

All I want to do is modify the existing kml file with:

 <href>http://developers.google.com/kml/documentation/images/top_left.jpg</href>
</Icon>
<overlayXY x="0" y="1" xunits="fraction" yunits="fraction"/>
<screenXY x="0" y="1" xunits="fraction" yunits="fraction"/>
<rotationXY x="0" y="0" xunits="fraction" yunits="fraction"/>
<size x="0" y="0" xunits="fraction" yunits="fraction"/>

How do I read and write into a kml file from a directory and make the modification from using this code? I am able to copy and paste the code into Google Earth  but  i want to be able to automate this process using python.


texts: Select polygon features that intersect a line feature using pyqgis

texts:
How do I select all features in a polygon layer that intersect (overlap) a line feature?

I need to implement this in python using QGIS, thus pyqgis. This code illustrates what I am looking for (line 10 is sort of pseudo-code, i.e. what I am looking for):

 1 area_layer_file = "./Areas.shp"
 2 line_layer_file = "./Lines.shp"
 3
 4 area_layer = QgsVectorLayer(area_layer_file, "Areas", "ogr")
 5 line_layer = QgsVectorLayer(line_layer_file, "Lines", "ogr")
 6
 7 areas = []
 8 for line_feature in line_layer.getFeatures():
 9     request = QgsFeatureRequest(QgsExpression(
10             area_feature.geometry().intersects(line_feature.geometry())
11         ))
12     areas = [area for area in area_layer.getFeatures(request)]

This image illustrates my question: the red line intersects the yellow selected areas.


texts: Google Earth Pro IOS app not opening exported kml file

texts:
I classified an image in QGIS using the dzetsaka plugin. I then exported it as a GeoTiff with the Rendered image Output mode selected, so that I could see it in Google Earth Pro. Otherwise, the raster only displays 0's and 1's. I made sure that my CRS was projected to EPSG: 4326 - WGS 84 since that is what Google Earth Pro is set to. I imported the file into GEP and used the Super Overlay process. Everything to this point works.

I then saved the file as a .kml and uploaded it into my drive to be uploaded into my GEP phone app (iphone 7). I am able to upload the file with the app, but it does not zoom into the correct position, nor shows what I see on the desktop GEP (image below). I tried creating a polygon within GEP on my desktop and my phone GEP app does recognize the file and works perfectly. I also noticed that my polygon kml is 2KB and my raster kml is 500 bytes.

I also tried vectorising my classified raster within QGIS, but because my resolution is so high, the process is taking 10+ hours and that is not sufficient (plus, it would crash on me every time). So I'd really like to stick with working with rasters.

Any suggestions?


texts:
I would like to pass multiple featureTypes to WFS protocol in order to accomplish a search feature in my app. the idea worked just fine with one featureType, but now I need to add another one to expand my search scope. Here is what I did:

var featureTypes = ["Layer_1", "Layer_2"];

var wfsProtocol = new OpenLayers.Protocol.WFS.v1_1_0({
             url: URL,
             featureType: featureTypes,
             featureNS: "http://mapmap.org",
             srsName: "EPSG:3857",
             version: "1.1.0",
             extractAttributes: true,
             isBaseLayer: false,
             visibility: true
               });

wfsProtocol.read({
    filter: filter,
    callback: processTheQuery,
    scope: [new OpenLayers.Strategy.Fixed({preload: true})],
         })

But still didn't work, Any idea?


texts: How to show "First Vertical Derivative of the Magnetic Field" survey data on a map?

texts:
I will be given a magnetic survey raster and I am being asked to display this on a map as a "First Vertical Derivative of the Magnetic Field". I've done this before with a different older dataset but I basically tried to match the raster colour ramp to a an example of a hard copy map.

I am not a geoscientist or geologist but I would like to understand (in a nutshell) what "First Vertical Derivative of the Magnetic Field" is and how to correctly symbolize it on a map.

I've seen some examples on the web and it seems fairly straight forward but, is there a predefined universal way of displaying this type of data?


texts: How can I accurately split roads from the TIGER roads dataset, to create proper junctions in a network dataset?

texts:
I am trying to create a network dataset in ArcMap using the TIGER roads dataset. A problem with this dataset is that some road segments are not split at intersections. Therefor, network analyst would treat it as though one could not turn onto that road from another (See picture below - the TIGER dataset is like the example of the "poor quality" in some cases)

I have used the "feature to line" tool to split roads where they intersect other roads for everything not classified as a primary or secondary road. However, I cannot do this for the primary or secondary roads, since it would likely be inaccurate to assume that you could come onto/off of the highway at any road that intersects. I have already tried using the Feature to Line tool for Primary and Secondary only where they intersect ramps, but it seems that in many cases one should be able to get onto the highway even where there is not segment classified as a ramp.

Has anyone use TIGER data to create a road network dataset before? Have you had this problem? Any suggestions as how to resolve it?


texts: Exporting cloud masked MODIS image collection from google earth engine

texts:
I am working on cloud masking, on a collection of Modis images in google earth engine.The code below works for the cloud masking part however when i try exporting the collection it returns an error as below
Error: Image.clipToBoundsAndScale, argument 'input': Invalid type. Expected: Image. Actual: ImageCollection.

I don't know where the problem is that is causing the error

var maskEmptyPixels = function(image) {
  var withObs = image.select('EVI').gt(0);
  return image.updateMask(withObs);
  };
var maskClouds = function(image) {
  var QA = image.select('SummaryQA');
  var bitMask = 1 << 10;
  return image.updateMask(QA.bitwiseAnd(bitMask).eq(0));
};
var table_bounds = function(image){
 return image.clip(table);
};
var collection =ee.ImageCollection('MODIS/006/MOD13Q1' )
               .filterDate('2019-03-01', '2019-07-30')
                .filterBounds(table).select('EVI','SummaryQA')
                .map(maskEmptyPixels);
var evicollection = collection.map(table_bounds);
var totalObsCount = evicollection
        .select('EVI')
        .count();
var collectionCloudMasked = evicollection.map(maskClouds);
var clearObsCount = collectionCloudMasked
        .select('EVI')
        .count()
        .unmask(0);
Map.setCenter(35.94,-0.37,8);
Map.addLayer(
    collectionCloudMasked.median(),
    {bands: ['EVI'],
     gain: 0.07,
     gamma: 1.4
    },
    'median of masked collection'
  );
  print(collectionCloudMasked);
var list = collectionCloudMasked.toList(11);
for(var i=1;i<11;i++){
  var nam=ee.String('EVI_img')
             .cat(ee.String(ee.Number(i)))
             .getInfo();
}
Export.image.toDrive({
       image: collectionCloudMasked,
       description: nam,
       scale:250,
       maxPixels:1e13
});


texts: Merge or dissolve line features based on two or more attributes

texts:
I am working in Arcmap 10.2 and I would like to do the following according to the example table below, but I do not have a clue on how to do this. I am working in ModelBuilder. Is there a tool for this or another easy way to do this?

The attribute table shows line features, some attributes have the same number (attribute B) and the same material (attribute A) this means that they are the extent of each other like the figure below:

If attributes have the same number and material I would like to make one Object of these ones and I would like it to sum the shape lengths. So for example ObjectID 1 and 2 will become one Object from Wood (Attribute A) with number 13324 (Attribute B) and a Shape length of 8.997 + 0.012 = 9.009. In the figure red and blue parts will be one part.


texts: Error while make Build osm2pgrouting on openSUSE 13.2

texts:
i am getting the error while make

/usr/lib64/gcc/x86_64-suse-linux/4.8/../../../../x86_64-suse-linux/bin/ld: cannot find -lpq
collect2: error: ld returned 1 exit status
make[2]: *** [osm2pgrouting] Error 1
make[1]: *** [CMakeFiles/osm2pgrouting.dir/all] Error 2
make: *** [all] Error 2

What i need to change?


texts: Selecting/tracing part of line feature using QGIS?

texts:
I have a shapefile which comprises a number of lines, each representing a bus route. I'm wanting to create a new bus route which is the same as an existing line to a certain point, but then branches off and follows a different route (see image below).

What's the best way to do this?

I'll be doing it numerous times. I'm thinking the best way is to draw a polygon around the parts I want and then clip using the same polygon, and dissolving the result.  I don't really want to draw new lines by snapping the vertices to existing routes as each route has thousands of vertices.

I am using QGIS 3.2.3.


texts:
I cannot add the below map server (Middle_East_Oil_and_Gas_Field_Map_by_CIA (WMTS) to qgis through add WMS/WMTS layer.

A connection will not be made to the following URL.

URL: https://tiles.arcgis.com/tiles/ULBqC49IEeIR01GF/arcgis/rest/services/Middle_East_Oil_and_Gas_Field_Map_by_CIA/MapServer?view=wmtsview&cacheKey=aec62228d3e2389d

What am I doing wrong?

QGIS version 2.18.26


texts:
I have a raster that I'd like to plot in the Azimuthal Equidistant projection with a grid and labels.  So far I've been able to do that.  However, I cannot get the labels inside the plotting area.  They are always cut off.

I have tried most of the par() options such as playing with par(mar), par(oma), and par(usr).  Hard to tell, but I think the par("usr") option does not work with raster::plot.  I also tried playing around with the extent of the raster, since R bases the plot size on the coordinate range of the initial data, but that also did not work.  Below is my code and cut off labels.

GSmean<-stack("[file][1]")
parallels <- c(80, 70, 60)
meridians <- seq(-180, 160, by = 30)
mlim <- c(-180, 180)
plim <- c(50, 80)
grat <- graticule(lons = meridians, lats = parallels, xlim = mlim, ylim = plim, proj = projection(GSmean))
labs <- graticule_labels(meridians, parallels, xline = 360, yline =46, proj = projection(GSmean))

#new plot
breaks<-c(-1,0,1,2,999)
cols <- c("darkgreen", "darkred","gray81","gray81")
plot(GSmean, col=cols,breaks=breaks,axes=F,legend=FALSE,box=FALSE)
plot(grat, add = TRUE, lty = 3,col="honeydew3")
text(labs, lab = parse(text= labs$lab),
     col= c("gray30", "gray30")[labs$islon + 1], cex = 0.85)


texts: Can QGIS open folder with "run feature action" tool?

texts:
Is possible to open folder with "run feature action" tool? I can open specific file like .dwg, but don`t know how open folder in specific location. This folder contain multiple .dwg files and when user click this action, folder open in another window.

Is this possible in QGIS 1.8.0, and how?


texts: Set missing pipe network data (age, size) based on nearest neighboring or connecting pipes

texts:
I have pipe network data (polylines) with missing age and size attributes.  For pipe age, I have 10779 pipes with 1034 of them having NODATA values.  For pipe size, I have 10779 pipes with 738 of them having NODATA values.

What tool can I use to set the missing values based on the nearest (or connecting) pipe attribute information?

EDIT 1: In my mind, I feel like I could set up something like this pseudo code:

for each pipe
  if the pipe[age]==NODATA
    find the pipe_nearest
    pipe[age] = pipe_nearest[age]
  if the pipe[size]==NODATA
    find the pipe_nearest
    pipe[size] = pipe_nearest[size]

And then it would be all done.  I just can't figure out how to determine the polyline nearest the current polyline.

Edit 2: Attempting to use FelixIP's solution.  First off, I got access to another copy of the data set with the pipes clearly connected to each other - I have from and to node fields (FNODE_1, and TNODE_1), and the pipes are clearly snapped to each other.  I have copied my initial layer twice, creating has_X and Missing_X layers for both age (YEAR_) and diameter (DIAM) information.

Starting with diameter first, I adjusted the script for the names of my own layers, so it reads:

def fillGap (d,fr):
 mxd = arcpy.mapping.MapDocument("CURRENT")
 lyr=arcpy.mapping.ListLayers(mxd,"Has_Diam")[0]
 q='"TNODE_1"=%s%s%s' %(r"'",fr,"'")
 tbl=arcpy.da.TableToNumPyArray(lyr,"DIAM",q)
 if len(tbl)==0:return d
 return tbl[0][0]

fillGap (!DIAM!,!FNODE_1!)

but for some reason I am getting this error on execution:

I am trying to sort it out, but since I don't completely the Python code that is used in calculating the expression, it is proving to be a little difficult.


texts:
I created a batch file, which executes a .fmw workspace on FME Desktop. I want this batch file to get started by the FME Server. So I created a simple workspace with a SystemCaller transformer, which calls the .bat.
Running that workspace on FME Desktop works perfectly.
As soon as I upload and run the workspace to the FME Server the logfile returns following error:

SystemCaller: Failed to Execute `\\*my UNC_Path*\start_fme.bat'

I granted a full access permission to the server and used a UNC-path to the .bat. On other projects I'm already using UNC-paths on FME Server workspaces to local data in the same directory as the .bat, so the connection between server and network drives is definitely working.
Any ideas, what I might have done wrong or might have forgotten?


texts: Cannot find polygons that are inside a big polygon using GeoPandas

texts:
First of all, some info about my system characteristics:

Python 2.7
GeoPandas 0.3.0

I am working on a program, where I filter out the data of two shapefiles -link with data- (each shapefile is a vector layer with polygons) and organize it according some statements. The two shapefiles look like the picture below, the picture is just a small part of the whole map.

After doing some filtering, I need to find in the remaining data those polygons inside the big polygons (in the picture the purple and the green polygon) and change some information if there is an intersection. I have tried using some of the GeoPandas functions for this task, but the results are not the ones that I expected.

Try 1:

# Read .shp files
b = gpd.read_file(fp + 'puchheim_buildings.shp')
la = gpd.read_file(fp + 'puchheim_landuse.shp')

#Small parts of the shape files
la_g = la.groupby('fclass').get_group('residential')
b_r = b[100:130]  # here I know there is an intersection with the big polygon

# Find the polygons inside big polygons
inter = gpd.overlay(b_r, la_g, how='intersection')

Here the work is well done, but the function returns me another shapefile and it takes a lot computing time if a use the hole data in the shapefiles.

Try 2:

# Read .shp files
b = gpd.read_file(fp + 'puchheim_buildings.shp')
la = gpd.read_file(fp + 'puchheim_landuse.shp')

# Find the polygons inside big polygons
inter = b.within(la.loc[0, 'geometry'])

In “try 2” I got a Boolean series with everything false, although I know that there are polygons inside the big ones. I tried to use instead the polygons their centroid and make the same calculations, but the result was the same, a Boolean series with all false

# Find the polygons inside big polygons
b_c = b.centroid
inter = b_c.within(la.loc[0, 'geometry'])

I would prefer a boolean series, because I just need to change some information in the shapefile if there is an intersection. The idea is something like this (the code below is a matter illustration from my idea):

# Read .shp files
b = gpd.read_file(fp + 'puchheim_buildings.shp')
la = gpd.read_file(fp + 'puchheim_landuse.shp')

# Find intersections
inter = b.within(la.loc[0, 'geometry'])

# Comparation
for r in range(len(b))
    if inter[r] == True:
        b[r]['type'] = 'residential'

How can I do that or have a better idea that achieves the result that I expect?


texts: Export raster with its min max range from GEE

texts:
I'm trying to export a raster as GeoTIFF from GEE and keep it's min max range (pixel value) so I can import the raster into ArcGIS and check it's min, max, mean and std.dev. in Symbology. My attempt so far can only manage to export the raster layer with a green color but when imported into ArcGIS all min, max, mean and std. dev. is 0.
Here is the code I have tried:

var Blueskin =ee.Image('COPERNICUS/S2/20160721T223713_20160722T000334_T59GMK');
var Red=Blueskin.select('B4');
var NIR=Blueskin.select('B8');
var ndvi = NIR.subtract(Red).divide(NIR.add(Red)).rename('NDVI');
var mean = ee.Number(stats.get('NDVI_mean'));
var std = ee.Number(stats.get('NDVI_stdDev'));
//find regions (Ulva) in ndvi has value higher than mean +1SD
var maskImage = ndvi.updateMask(ndvi.gt(mean.add(std)));
Map.addLayer(maskImage, {palette: 'green'}, 'mask image');
//get Ulva region max min
var stats = maskImage.reduceRegion({
  reducer: ee.Reducer.minMax(),
  bestEffort: true,
});

// Display min and max
print(stats);
//it's a bit confusing here but NDVI_min =maskImage_min according to print I just don't know how to change the name
var NDVI_min = ee.Number(stats.get('NDVI_min'));
var NDVI_max = ee.Number(stats.get('NDVI_max'));
print(NDVI_min);
print(NDVI_max);
var visualization_params =  {min: NDVI_min, max: NDVI_max, palette:'green'};
// Export a cloud-optimized GeoTIFF.
Export.image.toDrive({
  image: maskImage.visualize(visualization_params),
  description: '21072016UNDVI',
  scale: 10,
  region: table,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});

The min and max I set in the visualization_params doesn't seem to work and I can't use ee.Reducer.minMax() as an argument in visualization_params..


texts: Ubuntu Linux | Unable to configure ECW extenion in Geoserver 2.14.1 (war file)

texts:
We are using Ubuntu Linux 16.04, Geo-server war file 2.14.1 & Tomcat 8.5 (64 bit)

We have configured Geoserver in tomcat and trying to configure MrSID & ECW plugins

Steps Completed :-

1.GDAL ver. 1.9.2 installed

2.Geoserver Stable Plugins are added respective directory

3.Set an reference path's in Ubuntu Linux

4.Followed this steps from Geo-server site

5.Now MrSID extension works fine and for ECW extension shows error "terminate called after throwing an instance of 'std::length_error'
what():  basic_string::_S_create" & geoserver stops

6.Researching on this & it looks like its not supporting  for latest Linux version as pre this link (because if its not working for QGIS then doubt it will work for geoserver)

Tried this question and its steps

7.Is their any license required for this ? (reference site)


texts:
I am trying to crate a grid inside a shapefile, something like this. However, I am unable to generate such grid. I was wondering if anyone has an idea on how to accomplish this.

Here is my code -

WWWL.Shape<- readOGR("E:/Juan Arango", "WWL_Commerce_OK")
WWWL.Shape
plot(WWWL.Shape)
proj4string(WWWL.Shape)

bb <- bbox(WWWL.Shape)
cs <- c(3.28084, 3.28084)*6000  # cell size
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
grd
sp_grd <- SpatialGridDataFrame(grd,
                               data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(WWWL.Shape)))
plot(sp_grd)

Output of WWWL.Shape

class       : SpatialPolygonsDataFrame
features    : 1
extent      : 334367, 334498.7, 4088915, 4089057  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables   : 1
names       : Id
min values  :  0
max values  :  0

Output plot(sp_grid)

What I am hopping to get is something that looks like this


texts: Where can I find shapefiles for 1kmx1km grid (UTM grid) for the entire USA? It seems everywhere that I looked at there were shapefiles for states separately or for some particular UTM zones. I don't want to end up join 50 shapefiles.


texts: How to avoid overlapping SVG icons and labels?

texts:
I've dropped a point and then selected a SVG icon. But now I'm right next to the other points I want to label. Now you have the SVG icon, but be on the label, that is, the label may no longer be seen.

Is this a default setting that you can change?


texts: Modifying "pg_hba.conf" file to listen local TCP/IP conection?

texts:
I am starting with PostGIS and there is an error when trying to connect through a TCP/IP connection. What I have done is: installed PostGre 9.3 in my PC, then installed PostGIS 2.3 for postGre 9.3, then I tested the database in my computer using localhost and all worked fine, then I wanted to connect to the database form another computer that is in the same network (from now on, PC2).. I read many forums and then modify pg_hba.conf trying to allow PC2 to connect. I just couldn´t.

Default file:

# TYPE  DATABASE        USER            ADDRESS                 METHOD

# IPv4 local connections:
host    all             all             127.0.0.1/32            md5
# IPv6 local connections:
host    all             all             ::1/128                 md5
# Allow replication connections from localhost, by a user with the
# replication privilege.
#host    replication     postgres        127.0.0.1/32            md5
#host    replication     postgres        ::1/128                 md5

And this is how I modified it:

# TYPE  DATABASE        USER            ADDRESS                 METHOD

# IPv4 local connections:
host    all             all             127.0.0.1/32            md5
host    all             all             10.0.0.99               md5
# IPv6 local connections:
host    all             all             ::1/128                 md5
# Allow replication connections from localhost, by a user with the
# replication privilege.
#host    replication     postgres        127.0.0.1/32            md5
#host    replication     postgres        ::1/128                 md5

note that the new value is 10.0.0.99, which is the IP address of PC2.

Can anyone help?


texts: Sorting layers and attributes displayed in ArcMap Create Features window?

texts:
At ArcMap (advanced) 10.3.1, items displayed in the Create Features window are sorted alphabetically.

How can I re-sort them so that they match the Table of Contents display order?

A search through the Arc help documents did not discuss this issue...


texts:
I'm doing a project on vehicle tracking using GIS - for that I'm using the open source software package QuantumGIS (QGIS).

I would like to know how to customise QGIS using the Python console to edit the point shapefile geometries. I want the point file to get updated as the location of the vehicle changes.


texts: Reproject GEOTIFF tiles from UTM to Plate Caree

texts:
I need to reproject several GEOTIFF tiles(each having different height/width) in six different UTM zones from UTM to Plate Carree.

Output of gdalinfo of two sample tiles is enclosed. As is clearly seen they are the same type of projection but different UTM zone hence they are not the same projection.
 What I have been thinking of doing is reproject tiles from the other UTM zones to the tile of the first UTM zone and then mosaic all the tiles and then reproject to Plate Carree. I already have GIS software to do the reprojection and so my question relates to the validity of the approach.
Will this give the best result or are there other approaches that will give better results?

My ultimate goal is to cut tiles from the mosaic in such a way that height equals width.

Coordinate System is:
PROJCS["WGS 84 / UTM zone 42N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",69],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["meters",1],
    AUTHORITY["EPSG","32642"]]

and
2nd tile -

Coordinate System is:
PROJCS["WGS 84 / UTM zone 43N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["meters",1],
    AUTHORITY["EPSG","32643"]]


texts: Error validating geometry collection / multipolygon in PostGIS

texts:
I have just updated a table of polygons with a query which returned:

UPDATE parcel SET wkb_geometry = ST_Multi(wkb_geometry)

UPDATE 399212

Query returned successfully in 25 secs.

Now, all of the rows should contain multipolygons.

However, when I try validating geometries with ST_MakeValid, it turns out that there are geometry collections in stead of multipolygons.  How can I fix this or remove the geometry collections?

UPDATE parcel
SET wkb_geometry = ST_MakeValid(wkb_geometry)
WHERE ST_IsValid(wkb_geometry) = 'f'

ERROR:  Geometry type (GeometryCollection) does not match column type (MultiPolygon)
NOTICE:  Too few points in geometry component at or near point -121.78764087332523 37.700189714113364


texts: Azimuth calculation from One to Many Points in PostGIS

texts:
I have vector data of 3D City Model in PostgreSQL database.  So from my data geometry column, which is a surface geometry. (These are geometries of all the surfaces of buildings). I had extracted the points through ST_DumpPoints query as follows:

SELECT (ST_x(geom)) AS X, ST_y(geom) AS Y FROM (
     SELECT (St_DumpPoints(geometry)).geom from surface_geometry WHERE cityobject_id=95
   ) AS foo

Now I would like to calculate the azimuth from the origin (the below-created point) to all the points extracted with the above query from all the surface polygons.
 Considering this point created below as position, viewpoint or GPS position of an observer, I am building it with MakePoint as below. (So as observer position, I would be changing this point everytime or is there any way to automatise this?)

SELECT (ST_x(StandPoint)) AS X, ST_y(StandPoint) AS Y FROM (
     SELECT ST_SetSRID(ST_MakePoint(391139.27, 5821816.69), 25833) AS StandPoint) AS
 foo

I have seen in PostGIS that there is a way to calculate azimuth between two points. The function ST_Azimuth calculate between two points as in PostGIS manual.
But how to calculate azimuth from one to many?
Any tips on how to proceed?

Later I want to store that calculate the azimuth in a variable in python using PsyCopg 2 (This part is not difficult) for my further calculations.


texts: Setting coordinate system for image in ArcMap?

texts:
I am a beginner with ArcMap. I have a problem with setting coordinate system for my image. I am using an add-in which
creates shapefiles .shp (craters, area), grids, plots and etc. However I get errors (like the ones attached). Basically how do I
correctly set the coordinates so that ArcMap correctly interprets position of the shapefile on the map?

I know the number of columns and rows i.e. pixels, of my image and I know the resolution of my image (0.82m/px).
Link: http://wms.lroc.asu.edu/lroc/view_lroc/LRO-L-LROC-3-CDR-V1.0/M1142616950LC   (CDR PTIF (multi-resolution, lossy compression))

EDIT: I am attaching the tool which creates the shapefiles. In summary: with this tool one can measure the diameter and x,y of selected craters (circle shapes) and the area (polygons) in km^2 of the counting zone.


texts: Filter layer by pasting in WHERE clause?

texts:
I have a web map in Portal that has the Web AppBuilder Filter widget.

The map has a point layer in it called Work Orders.

I would like to filter the work orders by pasting in a WHERE clause.

The WHERE clause is generated in a separate work order management system.
Example: ((woclass = 'WORKORDER' or woclass = 'ACTIVITY')
and contains(description,'$sidewalk') > 0  and historyflag = 0 and
istask = 0 and siteid = 'SERVICES')

How can I do this?


texts: Using QGIS variables within SVG graphics

texts:
If you use SVG symbols within QGIS you can use params like these so you can adjust properties interactively in QGIS:

fill = param(fill)
fill-opacity = param(fill-opacity)
stroke = param(outline)
stroke-width = param(outline-width)
stroke-opacity = param(outline-opacity)

I would like to do the same with text within an SVG-image:
Use (text-)placeholders in the SVG and fill in the actual text stored as QGIS variables within the print composer:

Is this already possible?


texts: Seg Faults and Instantiation Errors while running PyQGIS on OS X

texts:
I'm having issues getting off the ground with running PyQGIS on my OS X. I'm frankly very confused as some errors occur sporadically while executing the same code.

When I try the following:

QgsApplication.setPrefixPath('/Applications/QGIS.app/Contents/MacOS', True)
QgsApplication.initQgis()
QgsProviderRegistry.instance().providerList()

I either get one of three scenarios:

Segmentation fault: 11 when returning the list
A blank list returned
Or a this error after initQgis()...

QCoreApplication::applicationDirPath: Please instantiate the QApplication object first
QWidget: Must construct a QApplication before a QPaintDevice

A workaround I found to work is to instantiate a QgsApplication object first as the error suggests with something like

app = QgsApplication([], True)

However, the problem here is that this creates a GUI which I don't need - I'm just trying to do some scripting. If I change the instantiation parameter to False, then I get an error with initQgis() of:

QWidget: Cannot create a QWidget when no GUI is being used

Also, in a lot of examples I've found online you don't need the instantiation step... is there something funky with my setup?

I'm using a venv if that makes any difference... all commands I want to run work normally in the QGIS Python console, it's just the external script that's giving me issues.

Any thoughts?

Python: 2.7.10,
QGIS: 2.10.1-Pisa,
OS X: 10.10.3


texts: KML file saved from KMZ missing tabular data

texts:
I saved a KML file in Google Earth from this KMZ ... https://www.rampp-team.com/documents/pennsylvania/watershed/LowerSusquehanna/lancaster_pa/Lancaster_CSLF.zip . Upon opening the KML file in QGIS, all of the tabular data is the same for all of the polygons, even though when the KMZ was loaded into Google Earth, there was color coding based on, apparently, tabular data.

Is there any way to preserve the data so that I can work with it in ArcView (shapefiles)?

There seems to be missing information in the KML that is present when displayed in Google Earth.


texts: Performing "mostly within" selection using ArcMap?

texts:
I'm looking for a way to determine if a polyline is "mostly within" a polygon.  Since mostly within is arbitrary, I would need be able to enter a percentage at its execution.

I thought there might have been an option in ArcMap using Select By Location, but there doesn't seem to be.  There are similar methods, such as "completely within" but I don't see how I can use any of them to achieve what I'm looking for.

Does this functionality or a method exist to achieve this?


texts:
I've got two tables which I want to combine and enrich with each others data. I'm getting stuck on the join part, as this needs to be done via a spatial query.

First the table structure: a table with addresses and a table with ensembles.

address: zipcode | number | addition | geom_point
ensemble: id | ensemble_id | geom_polygon

The resulting table I want to achieve has a row for every address with the corresponding ensemble id. Corresponding means: where the geom_point is within geom_polygon.

There can be addresses which have their point outside of the geom_polygon. The ensemble ID should then remain NULL

Is it possible to use some kind of spatial function in the join part of a query?

SELECT a.*, e.ensemble_id
FROM address as a
LEFT JOIN ensemble as e ON [[a.geom_point within e.geom_polygon]]


texts: Loss of Jupyter notebooks in ArcGIS Pro after installing Deep Learning modules

texts:
I've managed to break the ability to use Notebooks in ArcGIS Pro by installing the deep learning modules in a cloned environment.
These are the dependencies:
conda install -c esri -c fastai -c pytorch arcgis=1.8.2 scikit-image=0.15.0 pillow=6.2.2 libtiff=4.0.10 fastai=1.0.60 pytorch=1.4.0 torchvision=0.5.0 --no-pin

The cloned environment works fine in an IDE, but notebooks refuse to open in Pro.
I'm assuming one of the dependencies is the problem, is there a way to re-activate notebooks in Pro?


texts: Missing buffer values when using zonal statistics in ArcGIS Desktop?

texts:
I am using zonal statistics in ArcGis 10.2.2 to calculate the mean value of fire frequency within 1km buffers. However, I have noticed that this is not being calculated for all of my sites.  There are a few overlapping buffers in my dataset.

Is there a way to get around this?


texts:
I want to add a really simple legend to my mapview-map. I have a sfc_MULTILINESTRING-object called path and a sfc_POINT-object called places. I just wanted to have a legend that shows a line and says "path" and a point that says "Points of Interest". I dont really know how I could implement this in mapview. The bit of code at the moment looks like this:

mapview::mapview(path, lwd = 7) + mapview(places)


texts: ArcGIS API for Python Delete Multiple Features

texts:
Can you help me with my code below to delete/send multiple ObjectIDs to a delete parameter?

Currently the code just finds one record and send that Object ID to the delete paramter, but id like to be able to select a certain value in a attribute field and delete all or setup some sort of For or While Loop to pass all Object IDs to the delete paramter. One of the field names is 'Biologist' and a test record/value for the example is 'jah'.

from arcgis.gis import GIS
from IPython.display import display

gis = GIS(username="xxx", password="xx")

search_result = gis.content.search('test url')
search_result[1]

test_item = search_result[1]
test_layers = test_item.layers
test_layers

test_fset = test_layers[0].query()
test_fset.df

test_flayer = test_layers[0]
test_flayer.properties.capabilities

test_features = test_fset.features

remove_feature = [f for f in test_features if f.attributes['Biologist'] == 'jah'][0]
dele = remove_feature.get_value('OBJECTID')

delete_result = test_flayer.edit_features(deletes=str(dele))
delete_result


texts:
I have a quick editing question.
 I have a main line (street) from which I want to extend other small lines to reach some points (manholes), but I want to keep just one feature (the main line). I can't use the extent tool, because it works only with line-line case.

ps: I am using ArcGis 10 (unfortunately;P )

Thanks,
Ligia


texts:
Similar ESRI arcgis.com can anyone give pointers on how to move layer up/down? And if possible include the context menu.

Using JavaScript API


texts: How to calculate the average distance among set of points as measure of closeness

texts:
I have three series (or tables) of points in PostGIS which I want to compare. I would like to know how "close" each series is when compared to the other two. In other words I would need some measure such as the average of the average distance of each point with all the others (that is, the average of the distance matrix)... Sounds complicated, but I am sure some common statistic exist to compare series of points.

How my query would look like?


texts: What is the fastest way to iterate through a list of features?

texts:
I have written a QGIS plugin and it iterates through all features in a layer using provider.select() and compares geometry.  The problem, the plugin is very slow.  If there are 114 features in a shapefile and it needs to iterate 80 times then it will have to repeat the code 9120 times.

How can I speed up the iteration process?  Should I store all features in a list and then iterate thru the list?  Store all features in a memory layer (not sure what that is yet)?

Thanks

Ed


texts: Converting large amount of points to polygons with QGIS

texts:
I have a large point dataset of about 4 million points in QGIS. I would like to convert these to polygons by connecting the outer edges and hence create ribbon-like X shapes (see image below). I checked some tools and the concave hull is exactly what I am looking for. However, QGIS gets stuck because of the large extent of the dataset.

I tried to give all points a cluster ID with DBSCAN and do the concave hull per cluster ID, but the amount of clusters is too high to do this manually (about 18,200).

Is there another tool that can do the same but requires less computational efforts?

I took a small subset of the data and the result should look like this:


texts: Copying value from file geodatabase field into shapefile field using ModelBuilder?

texts:
Does anyone know how I can pass the value from a field in a file geodatabase to a field in a shapefile using ModelBuilder?

I have a Model to do a Definition query via the Select, then use Summary Statistics to SUM a column and write it to a file geodatabase table. I now need to take this table field value and write it to a shapefile field.