texts: How to handle large tables in QGIS using python?

texts:
I'm currently writing a plugin that calculates values based on a QgsVectorLayer's attributes. To display and print the values I'm using a QgsComposerTextTableV2. I have several problems with my approach:

QGIS crashes when I try to run my plugin (virtual function call or instant minidump - I think it has to do with QgsComposerView)
I need to display/print 60+ pages, but after 10+ QGIS nearly stops responding
if QGIS doesn't crash and I try to create a PDF, it takes more than an hour for ~60 pages
I couldn't find a way to increase the height of the table headers
the value of one cell can consist of multiple lines, that's why I'm leaving the frame calculation up to QGIS

Is using QgsComposerTextTableV2 a good idea or are there better options to display/print large tables in QGIS/python?

edit: ReportLab seems to be a good option, testing it now.

QGIS version 2.14.3

Below is a simplified version of my code, which can be executed in the QGIS python console

from qgis.core import QgsComposition, QgsComposerTextTableV2, QgsComposerFrame, QgsComposerTableColumn
from qgis.gui import QgsComposerView
import qgis.utils

windowTitle = 'myComposer'
iface.createNewComposer(windowTitle)  # returns virtual that can't be used

# composition
mapSettings = iface.mapCanvas().mapSettings()
composition = QgsComposition(mapSettings)
paperHeight = 297
paperWidth = 420  # 420 results in 422 somehow using QPrinter, QPainter directly with QComposition
composition.setPaperSize(paperWidth, paperHeight)
composition.setPlotStyle(QgsComposition.Print)

# table
textTable = QgsComposerTextTableV2(composition, False)
textTable.setWrapBehaviour(1)
textTable.setHeaderMode(1)  # show header on all pages
textTable.setResizeMode(3)  # RepeatUntilFinished
composition.addMultiFrame(textTable)

# frame
x = 20
y = 10
frameHeight = paperHeight - 2*y
frameWidth = paperWidth - x - y
tableFrame = QgsComposerFrame(composition, textTable, x, y, frameWidth, frameHeight)
textTable.addFrame(tableFrame)

# columns
columns = []
x = 10
i = 0
while i < x:
    tempColumn = QgsComposerTableColumn()
    tempColumn.setHeading('heading that doesn\'t fit\n into a single line')
    tempColumn.setWidth(35)
    columns.append(tempColumn)
    i += 1
textTable.setColumns(columns)

# content
contentList = []
x = 10
y = 500
i = 0
j = 0
while i < y:
    newRow = []
    while j < x:
        newRow.append('row ' +  str(i) + ' column ' + str(j))
        j += 1
    contentList.append(newRow)
    j = 0
    i += 1
textTable.setContents(contentList)
textTable.recalculateTableSize()

# get composerView
activeComposers = iface.activeComposers()
for composer in activeComposers:
    if composer.composerWindow().windowTitle() == windowTitle:
        composer.setComposition(composition)
        break


texts: Mapserver+openlayers: displaying a shapefile with Cylindrical Equal Area (Lambert) projection

texts:
I have some shapefiles in "World Cylindrical Equal Area (Lambert) projection".
The .prj reads:

PROJCS["World_Cylindrical_Equal_Area",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_
1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PR
OJECTION["Cylindrical_Equal_Area"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",
0.0],PARAMETER["Central_Meridian",-160.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1
.0]]

In consequence, I wrote the following mapfile:

MAP
    NAME "mapTest"
    OUTPUTFORMAT
      NAME "png"
      DRIVER AGG/PNG
      MIMETYPE "image/png"
      IMAGEMODE RGBA
      EXTENSION "png"
      FORMATOPTION "GAMMA=0.75"
      TRANSPARENT ON
    END
    SIZE 800 600
    EXTENT -20037507.067162 -6363885.331926 20037507.067162 6324287.725856
    UNITS meters
    SHAPEPATH "/data/"
    IMAGECOLOR 255 255 255
    FONTSET "/usr/share/fonts/truetype/ubuntu-font-family/bruno_fontsets.list"

    WEB
      IMAGEPATH "/data/tmp/"
      IMAGEURL "/tmp/"
      METADATA
        wms_title "WMS Demo Server"
            wms_onlineresource  "http://localhost/cgi-bin/mapserv?map=/var/www/ex/mapTest.map&"
        wms_enable_request "*"
      END # metadata
    END # web

    PROJECTION
      "proj=cea"
       "lon_0=-160"
       "lat_ts=0"
       "x_0=0"
       "y_0=0"
       "ellps=WGS84"
       "datum=WGS84"
       "units=m"
       "no_defs"
    END

    #######################"
    LAYER
      NAME "countries"
      STATUS ON
      TYPE POLYGON

      PROJECTION
        "proj=cea"
        "lon_0=-160"
        "lat_ts=0"
        "x_0=0"
        "y_0=0"
        "ellps=WGS84"
        "datum=WGS84"
        "units=m"
        "no_defs"
      END

      DATA "/data/data_repo/WRI/Base_Data/Countries/Countries"
      LABELITEM "COUNTRY"
      PROCESSING "LABEL_NO_CLIP=ON" #avoid repeating labels
      CLASS
        NAME "countries"
        STYLE
          OUTLINECOLOR 200 200 200
          WIDTH 0.5
          ANTIALIAS true
        END # style

        LABEL
          COLOR 50 50 49
          MINFEATURESIZE 20
          MINDISTANCE 10000
          REPEATDISTANCE 0
          BUFFER 5
          SIZE 12
          FONT "Ubuntu-C"
          TYPE truetype
          POSITION auto
          PARTIALS FALSE
        END # label

      END # class
      METADATA
        wms_name "countries"
      END # metadata
    END # LAYER
    ###########################

END # MAP

and the following openlayers code:

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">

Unfortunately, the displayed image is completely messed-up:

Can someone find what is incorrect in my codes?


texts: Creating a vector layer that fills the space between two line vector

texts:
I need to create a polygon vector that fills the space between two line vectors I have.

The black line vector represents the coastline, while the blue line vector indicates the 250m depth contour.

The idea is to create a polygon that represents the entire marine territory of the region that is under 250m depth (meaning, between the two line vectors)


texts: Converting LAS point cloud from feet to meters in R

texts:
Building off of a previous question:
Clean way to convert from US feet to meters with lidR - LiDAR point cloud
I am working with the USGS LiDAR point cloud data and am still having this rescaling or whatever you want to call it using the:
  las$X = las$X * ft2m
  las$Y = las$Y * ft2m
  las$Z = las$Z * ft2m

Approach from above. I get the error:
Error: Trying to store values ranging in [3496056, 3497580] but storable range is [7852516.353, 12147483.647]

My Z and Y coordinates are no issue, it is only the X. What's up, because it is driving me insane. Only four of my tiles seem to have this problem.


texts: How to tranform MODIS tiles into lat/long?

texts:
I'm trying to get the vegetation index nearest to a particular lat/long. I'm getting the data from here which is 250-meter spatial resolution hdf file. The globally gridded projection is split in 36 (horizontal) by 18 (vertical) tiles, each with a resolution of approximately 10° by 10°.

Now each tile has 4800x4800 values. I can extract the GRINGPOINT(GRINGPOINTLATITUDE/GRINGPOINTLONGITUDE) from the metadata but I'm not sure how to store values of the VI(vegetation index) corresponding to the lat/long inside the tile.

Let us suppose the GRINGPOINT is like this:

GRINGPOINTLATITUDE = [-0.0068357003079556, 9.99897831672069, 9.9909309627606, 5.67994760615426e-06]
GRINGPOINTLONGITUDE = [-179.999951582871, 179.928473473918, -169.920147289013, -169.99173290556]

Which value should I take in place of latitude and longitude in the below code ?

latitude = -90.0
longitude = -180.0
count = 0
res = []
for i in range(0,4800):
    del res[:]
    for j in range(0,4800):
        data = {"lat":latitude,"lng":longitude,"ndvi":values[i,j]}
        res.append(data)
        longitude += 0.00225
    collection.insert(res)
    latitude +=0.00225
    longitude = -180

If there is any other way of doing this then please let me know.


texts: Spatially Joining Many Shapefiles. Use Spatial Join or UpdateCursor?

texts:
I have a large number of tree records in 100s of point shapefiles. Each shapefile contains a number of trees that fall within 10-50 parcels. The parcel layer contains ~60,000 shapes. For each point shapefile, I need to extract some of the features from the parcel layer and add them to the point layer.

I have two workflows in mind, none of which I'm very excited about and would like feedback. Working in ArcMap 10.2 advanced.

Spatial Join:

Do a spatial join for each tree layer with the parcel layer. Since the output will have to be further processed (more fields will be added, and then output to an excel table), I don't like that spatial join creates a new shapefile, and adds fields that I don't need in the final output. I have also found spatial join to be rather clunky.

Using a nested search cursor and update cursor. (See Using Select By Location to update field in feature class using ArcPy?).

Basically, da.SearchCursor(Parcels), SelectLayerByAttribute(Parcels), SelectLayerByLocation(Trees), da.UpdateCursor(Trees), create new fields and populate them with values from the selected parcels.
The problem here seems to be that I'm looping through the entire parcel layer despite knowing that only a few parcels need to be considered.

The use of an in_memory table of only the parcels I needed was what I was missing.

Here's the code I ended up using:

import arcpy, os

#Set environment and define variables
arcpy.env.workspace = r'C:\...\Example'

path = r'C:\...\Example'
outpath = r'C:\...\Example\SurveyTrees.gdb'
parcels = r'C:\...\Example\ParcelProcess.gdb\Parcel'

#Make a list of all the shapefiles in directory
shapes = [os.path.join(path,shp) for shp in arcpy.ListFeatureClasses('*')]
outname = os.path.abspath(shapes[0]).split(os.sep)[-3]

#Merge shapefiles
arcpy.Merge_management(shapes, os.path.join(outpath, outname))

trees = os.path.join(outpath, outname)

#Make feature layers
arcpy.MakeFeatureLayer_management(trees, "tree_lyr")
arcpy.MakeFeatureLayer_management(parcels, "parcel_lyr")

#Select parcels that intersect trees
arcpy.SelectLayerByLocation_management("parcel_lyr", "INTERSECT", "tree_lyr",
                                   "#", "NEW_SELECTION")

#Create in memory feature class of selected parcels
arcpy.FeatureClassToFeatureClass_conversion("parcel_lyr", "in_memory", "treeparcels")
treeparcels = r'in_memory\treeparcels'

#Make feature layer of relevant parcels
arcpy.MakeFeatureLayer_management(treeparcels, "treepar_lyr")

#Create search cursor for relevant parcels
with arcpy.da.SearchCursor("treepar_lyr",
                       ["OID","LocName", "SITE_ADDR",
                        "CITY", "ZIP"])as pcur:
    for prow in pcur:
        #Select one parcel at a time
        arcpy.SelectLayerByAttribute_management("treepar_lyr", "NEW_SELECTION",
                                            "OID = {}".format(prow[0]))
        #Select trees that are within that parcel
        arcpy.SelectLayerByLocation_management("tree_lyr",
                                           "WITHIN", "treepar_lyr",
                                           "#", "NEW_SELECTION")
        #Create update cursor for trees
        with arcpy.da.UpdateCursor("tree_lyr",
                               ["Loc_Name", "SiteAddr",
                                "City_", "Zip_Code"]) as tcur:
            #Update tree attributes with attributes from selected parcel
            for trow in tcur:
                trow[0] = prow[1]
                trow[1] = prow[2]
                trow[2] = prow[3]
                trow[3] = prow[4]
                tcur.updateRow(trow)


texts:
I'm working with Leaflet and I'd like to show the day and night zones adding a dynamic layer to my map. I was thinking to add a static layer with some transparency and make it move every x minutes but it's tricky because the shape of the shadow... changes all the time during the year.

Does anyone know if there exists some open source layers for every season or better if I could find a way to calculate an accurate geojson layer for the earth's shadow ?


texts: Method to join line segments in a road network that are shorter than X length to an adjacent segment using ArcGIS?

texts:
I'm using ArcGIS Desktop 10.2 and I'd like to take roads in a network shorter than a certain length and append them to an adjacent street section. Currently the road network files are split intersection to intersection which creates awkwardly small segments.

Could you suggest ways of automating the process to select all roads under a certain length and append them to an adjacent street section?


texts: Manipulating data with a Handsontable spreadsheet in a Leaflet GeoJSON Layer and rerender the map

texts:
The goal is to add a spreadsheet view to a GeoJSON Layer to comfty edit the properties directly in the spreadsheet and the location (lat/long) within the map by dragging the marker with leaflet.draw.
For the speadsheet functionality i found a nice lib called Handsontable:
https://github.com/handsontable/handsontable
Handsontable works with array so I push the properties data  with onEachFeature in the array tableEditArray and in the Leaflet Layer called drawnItems and bind a mouseover event on the marker to view the data:
onEachFeature: function(feature, layer) {
                       tableEditArray.push(feature.properties);
                       drawnItems.addLayer(layer);
                       layer.on('mouseover', function_to_view_properties);
                }

Then I add drawnItems to the map and initialise a Handsontable with the data from the array tableEditArray.
So far, so good, the map is drawn and the spreadsheet is editable:

If the user edits a cell in the sheet, how can I update the drawnItems Layer and rerender the map?


texts: con statement to change the default value in ArcGIS

texts: I have obtained a flow accumulation map using Arc Hydro tools and obtained the default value in the range of 0-204567. Now I need to change these value by default so that it gives 1 instead of zero. I have tried to reclassify it but I am not getting the range 1- 204567. So for this I have to use a con operation in raster calculator but I don't know what condition to use so that I get a value of 1-204567. How do I specifiy the correct con statement in order to change the value?


texts:
If I am going to use cellstat alone, CellStat would produce mean of 2000 rasters inside the FOLDER(ALL RASTERS). I want to average only every 5 raster which would result to 400 rasters.

How do you do this in ArcGIS 10.1?


texts: How CARTO (or other GIS-software) decides where labels are located

texts:
I´m calculating the distances between different places in my region/country.

I don´t know where I can find a database of all the towns/cities already geolocated, so I have a question for considering a "middle point" of each region of the shapefile:

How CARTO (or other GIS-software) decides where to put the labels (in bold in the image below) of each feature?

My target is getting the centroid and consider that that´s the center of each town/city.


texts: Zonal Statistics as Table gives error when there is no overlap

texts:
I'm trying to generate Zonal Statistics with a raster and set of polygons.

In my case the polygons represent census tracts in the US and the raster is population density.

The problem is that sometimes the raster and polygon don't overlap, because the raster only includes the lower 48 states, and I'm writing a script to loop over files from all states to do this operation, and those files include some Alaska/Hawaii files.

Is there a function to check if a raster and a polygon overlap?

Right now I'm working on ArcMap and just copying as snippet as I go:

This gives an error, because state number 02 is Alaska. Ideally it would be smart enough to skip this case.


texts: Extract highest quality MODIS pixels in Earth Engine

texts:
Based on the MYD09GA QC_500m layer description below:
MYD09GA LAYERS info QC_500m: 500m Reflectance Band Quality
Bit 2-5: band 1 data quality four bit range
0000 (0): highest quality

I need to extract only the highest quality data of band1. My current code is given below as well. After applying the mask, when I try to generate a chart of the timeseries for Band1, I receive the error:
Error generating chart: Data column(s) for axis #0 cannot be of type string.
Is something wrong with my mask function? or I need to do any data transformation for the timeseries generation? Can anyone help out please. Many thanks

var mod = ee.ImageCollection("MODIS/006/MYD09GA"),
geometry1 = /* color: #d63000 */ee.Geometry.Point([10.644168108701706, 64.4188479841206]);

//Selecting Band 1 and QC layers
var dataset = mod.filterDate('2013-01-01', '2013-12-31').select('sur_refl_b01','QC_500m')
var single = dataset.first()

//create mask to extract only highest quality band 1 data
var filter = function(image){
var mask = single.select('QC_500m').eq(0)
return ee.Image(image).updateMask(mask)
}
var maskapply = dataset.map(filter)
var chart = ui.Chart.image.seriesByRegion({
imageCollection: maskapply,
regions: geometry1,
reducer: ee.Reducer.mean(),
scale: 1000,
band: 'sur_refl_b01',
})
print(chart)


texts: How can I get the original DEM from a hillshade raster (GeoTiff)?

texts: The problem is that the hillshade raster has a zfactor of 1.4 and for hydrological tool I need the real data. I tryed to divide the whole raster with raster-calculator with a factor of 1.4 but it cannot creat the same original data. Is there any tool for doing that?


texts: Summarizing point data to multi point data and creating convex hull in R

texts:
I have latitude and longitude information like 26.847905 73.9575005 for about 800 points, I need to trim these lat/lon information to say 3 decimal places.

I used the following syntax in R

# Rounding off Coordinates of each of point to 3 decimal places
shp1_DE$lon <- round(shp1_DE$coords.x1, 3)
shp1_DE$lat <- round(shp1_DE$coords.x2, 3)

Now I have many points with lat/lon for example 26.847, 73.957 etc...

I want to summarize data to these points and form point-multipoint shapefile.

After that I need to create convex hull for each of these features. How to go about it in R?


texts: Why do first returns have more blanks than last returns in LiDAR data?

texts:
I'm working on an ASCII LiDAR file that has the first and last pulse returns included.

Can anyone tell me why a first pulse return may come back with more blank spaces in the data than a last pulse return? I just think logically it should be the other way around, but maybe I'm missing something?


texts:
I work actually on Network Detection. Precisely, I am detecting public light electric cables undergound and drawing a map with lines connected on GPS points. Each of these points gets XYZ, and the depth of the cable, according to surface altitude.

For quality control, I have to check if there isn't more than 35 cm between to points, because it means that the detection accuracy isn't good enough.

My complete plan on QGIS look like this (fictional exemple) :

Currently, i check visually if there isn't a 35 cm depth gap between two points. I would like to make a model builder that could check automatically.

I need to make QGIS understand that he must calculate (A depth - B depth) then (B depth - C depth) ect ect... in that specific order.
It implies to make it understand he has to follow the points along the line in the start/end order.

I cannot follow XY because here there is only one line but it's an exception. Generally, there's a lot of lines that cuts or follow themselves so it's a big mess.

I ran out of idea. I checked others topics about a similar subjects but couldn't find an answer.

Thanks.


texts: QGIS - import *.dsf files (Danish National Format)

texts:
I'm new to QGIS and in way above my head with the quest at hand. I'm trying to create a 3D map of a Danish city, to this task the Danish municipality has provided a map of the area the so called KORT10 which I would like to import into QGIS, the map consist of four parts which all is DSF files, it should be vector format, but when trying to add vector layer it pops up with an error saying:

"Invalid Data Source C:\Users\.........\KORT10.dsf is not a valid or recognized data source".

I have heard some talk about it being Shapefiles or something like that, but as mentioned this is the very first time using QGIS and GIS in general.

I have uploaded the file from the municipality if what I have just described made absolutely no sense but somebody has an interest in helping me. The CRS system used is "EPSG2583/ETRS89"

http://www.filedropper.com/510dsflutm32-euref891

I hope somebody has the time and knowledge to help me.


texts: Sum Line length and repeat lines in QGIS?

texts:
I have a grid of 10 km x 10 km polygons, and am planning survey transect lines across the area.  I am trying to use the Sum Line Lengths tool to calculate the total transect length in each polygon.  The transect line is a single, continuous line, which zig-zags across the grid.  The CRS is UTM, units metres.

The issue I have is that the summed length is incorrect - it is less than it should be.

I think this is due to the transect line repeating its track across the grid cells - might this be the case?

So the transect crosses over precisely the same ground a few times.

Another minor issue, but perhaps a clue to what is happening, is the line count is also incorrect: it's either 0 or 1, when in some cells it should be 2 or 3.

I've done various tests with line layers on the same grid of polygons, and Sum Line Length gives the correct output if the lines, or single line with multiple passes through the cell, do(es) not repeat its track precisely as I want it to in reality...  The number of lines / cell is always wrong though.

I am using QGIS v2.18.16.


texts: Identify area within one shapefile not in another shapefile in R

texts:
I have two shapefiles. I would like to identify the part of the first shapefile that doesn't overlap with the second shapefile. Then I would like to clip that part of the first shapefile and make it its own shapefile to merge with the second shapefile in R.  I have tried intersect and union in the raster package and gIntersection in the rgeos package but haven't had any luck.

After trying to use the answer below from Simbamangu I think the reason the answer is not working for my shapefiles is that the coordinates for my two shapefiles are formatted differently.  The answer works with the example data.

Shapefile attributes

Shapefile 1 output from attributes             Shapefile 2 output from attributes

  $bbox                                            $bbox
        min       max                                min       max
x -67.33333 -66.33739                              x  639772.8  723467.9
y  40.58634  41.50000                              y 4499609.0 4597064.1

Is there a way to format the coordinates so both shapefiles are formatted the same way and could this be the reason gIntersection is not working.

intersect1<-rgeos::gIntersection(shapefile1, shpaefile2)
NULL

I figured out to format the coordinates from the 2 shapefiles into a similar format.  The gDifference function doesn't seem to be identifying all of the differences between the two shapefiles.  See output below.

Output from getting answer to work.  The only part of the Shapefile 1 identified as outside of Shapefile 2 is the little sliver of green, but the 2 sections in white should also be identified.

Here are two example polygons - not the actually shapefiles

# Create SpatialPolygons objects
polygon1 <- readWKT("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20))")
polygon2 <- readWKT("POLYGON((-180 -20, -140 70, 10 0, -140 -60, -180 -20))")

R information

R version 3.3.1 (2016-06-21)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] grid      datasets  utils     stats     graphics  grDevices methods
[8] base

other attached packages:
 [1] taRifx.geo_1.0.7    devtools_1.12.0     rgeos_0.3-19
 [4] raster_2.5-8        calibrate_1.7.2     spatialEco_0.1-5
 [7] spatstat_1.46-1     rpart_4.1-10        nlme_3.1-128
[10] rgdal_1.1-10        maptools_0.8-39     mapdata_2.2-6
[13] maps_3.1.0          sp_1.2-3            car_2.1-2
[16] xlsx_0.5.7          xlsxjars_0.6.1      rJava_0.9-8
[19] plyr_1.8.4          MASS_7.3-45         RODBC_1.3-13
[22] latticeExtra_0.6-28 RColorBrewer_1.1-2  lattice_0.20-33


texts: Scattering overlapping points using ArcGIS Desktop?

texts:
I am working with a series of points that are located at the same location.

Can I scatter these point so that it is possible to visualize them as I zoom in?

Unfortunately, they don't have any element in common in the Attribute Table.

I am using ArcGIS Desktop.


texts: Adjusting text border thickness in MapInfo

texts: MapInfo allows you to add an outline to your text by using the 'border' or 'halo' function. Is there any way to customise the thickness of these outlines? I find the standard setting far too thick, and adds to much clutter to my maps.


texts: Split a polylines according to a specific point in R

texts:
I have 2 distinct spatial objects:

on one hand I have a SpatialLinesDataFrame
on the other hand I have a SpatialPointsDataFrame

Considering one line and one point that intersects that line, I would like to split the line into two segments according to the point.

The point is fixed and does intersect the line, the following command from the rgeos package

gIntersects(line.shp,point.shp)

returns TRUE.

A similar question has been asked here but the answer is based on the st_split function. However, this function is no longer implemented in the latest version of the sf package (version 0.7-3).


texts: how Npgsql.EntityFramework support Spatial data

texts:
I am using Npgsql.EntityFramework + EntityFramework6.My database is postgresql. I Added to my C# project with this tutorial :

  https://github.com/npgsql/npgsql/wiki/Visual-Studio-Design-Time-Support---DDEX-Provider

It works well. The problem is this that when I have a field of spatial data like point and polygon it cannot add that field to model. How can I solve this problem?


texts: Geoserver WMS Filtering: Clipping geometry against a polygon

texts:
My data set is a series of multilinestrings representing roads. During a WMS GetMap request, I want to specify a polygon and display only the portions of those roads that fall within the polygon. I am fine using ECQL or SLD.

My current CQL_FILTER looks like:

INTERSECTS(wkb_geometry, POLYGON((-117 34 ... ...)))

This returns all roads intersecting the polygon, but displays the full extent of those roads - i.e. it displays all vertices/segments for the entire road, not just those contained within the polygon. Image example:

Is there a way to extract the vertices from each geometry, and filter against those individual vertices, instead of the entire line?

Something like the following command, although this doesn't work,

INTERSECTS(VERTICES(wkb_geometry), POLYGON((... ...))))

Failing with the error message:

  Rendering process failed Binary geometry filter, but first expression is not a property name? (it's a class org.geotools.filter.function.FilterFunction_vertices)

This makes sense, as VERTICES(wkb_geometry) doesn't actually return a property of the dataset, but an entirely new set of data. I assume this is called "geometry transformation", but I can't figure out whether it even CAN be done in this way, and if so, how.


texts: Openlayers-3, can't get attributes from a WMS layer

texts:
I have an osm map and a WMS from mapserver in Openlayers-3. I am trying to display some attributes from WMS "emprises" layer according to this example. But The window is empty : I just see "Feature1234:'' ".
I have add PropertyName in the event map.on, but dont'know how to call/write it in the iframe part...if I am in the right direction !
WMS request is enabled (template) in the mapfile and works fine with Qgis.

<link rel="stylesheet" href="http://openlayers.org/en/v3.17.1/css/ol.css" type="text/css">

&nbsp;


texts: Visualize Spatial Data in Reporting Service

texts:
Is it possible to visualize a piece of geographical map from SQL server Spatial data (geometry, not geography) in a report from reporting service? If yes, how should it be done in a course of action.

A concrete example is to display the orange line from the Chinese subway's map in a specific reporting page.


texts: How to choose an area (of mxn metres) within a raster that have values above a certain range

texts:
I've got a raster of concentrations of iron on Mars and I want to be able to choose several areas of either m by n dimensions in metres, or just any area in m^2, that encapsulate regions that have high concentrations.

For example:  I want all regions with concentrations over 5 for an area of 50x50 metres, or an area of 2500 m^2.

What process would you use?  Would I create a grid with 50x50 metre cells? Then perform some sort of calculation that tabulates the average cell concentration and output concentrations over some value?  I don't have any idea how to do that, but even if I did the areas don't need to be perfect squares so it wouldn't be the  best solution just better than not do any analysis.  Either way if anyone has any ideas how to choose areas like this either way I'd really appreciate any help.

Thanks


texts: Getting shapefile in right projection using ArcGIS Desktop?

texts:
I have an older map that I didn't create. I'm just trying to get everything to line up on it, and having a hard time.

The data frame is in NAD_1927_UTM_Zone_12N, as are the point shapefiles (derived from a .dbf table that I had to repair the data source to) and one line shapefile.
All the other data is in StatePlane (either 1927 or 1983) and that's what shows up in the layout view (aerial, roads, buildings, polygons, etc.)

I've changed the points (via right click, data, repair data source, browsing for the data and then editing "Coordinate System of Input Coordinates" to StatePlane, but that doesn't work. The point shapefiles aren't visible anywhere, even if I zoom to the layer or go to Full Extent.

When the points are in the UTM projection, they're at a ridiculous scale (1:1,000,000) and they're all in a straight line instead of spread out like the old PDF version of the map I have.

I'm trying to make some simple edits, but I can't fix the bigger issue of shapefiles not lining up.

Although I don't know if this is the problem, when I open the map I get a notice that says, "This document may contain VBA code but you don't have a VBA license." I would just assume that scripts don't run, not that that would affect projections...

Any insight would be greatly appreciated. Thanks.

edit: No, the pdf is uneditable. I did realize the XY were being mapped to the same field, which was wrong (the XY data is in northings and eastings), but correcting that didn't fix the problem. Unfortunately, whenever I try to add data, browse for projections and the like, ArcMap crashes which is an unrelated server problem. I could copy the data to my own folder but there are a lot of files and some of them are linked and from other folders and ...still of course in different projections. Although I fixed the VBA license problem, so that has nothing to do with it.

I tried doing Data Management > Projections & Transformations > Features > Project and each time (regardless of which projection I choose) the new points project exactly where they currently are. I'm currently trying suggestion #1 at the moment.


texts: Why does AddGeometryColumn error using the libspatialite-2 DLL with pysqlite2?

texts:
What I'm trying to achieve is to create a geometry column of a point based on the latitude and longitude of that particular row in the already created table.

Given the following code:

#!/usr/local/bin/python
from pysqlite2 import dbapi2 as sqlite3

conn = sqlite3.connect('adde.sqlite')
conn.enable_load_extension(True)
conn.execute("SELECT load_extension('libspatialite-2.dll');")
conn.execute("SELECT InitSpatialMetadata();")
conn.enable_load_extension(False)
conn.execute("SELECT AddGeometryColumn('coords','geom_column', 4326, 'Point', 2);")
conn.execute("UPDATE coords SET geom_column = GEOMETRYFROMTEXT('POINT('||longitude||' '||latitude||')', 4326);")
conn.close()

I receive the error:

AddGeometryColumn() error: "insert on table 'geometry_columns' violates constraint: 'spatial_ref_sys.srid'"

I've been able to create the column in the Spatialite GUI no problem, but I'd like to do all of this programmatically using Python.  I've searched for a while now and I'm coming up a bit short.  Any guidance as to why I'm receiving this particular error would be greatly appreciated.

I'm using Python 2.6 and pysqlite2 on Windows 7 64-bit


texts: How to display geotagged photos in layout view to export in map book ArcGIS for Desktop?

texts:
Im trying to create a 140 page map book using data driven pages and the production mapping extension in ArcGIS 10.2.2 and/or 10.3. As the image shows, the production mapping extension has a Graphic Table Element which allows me to display the chosen attributes of visible features on each page (driven by a grid index). I have a feature class of geo tagged images with path attributes.

I would like to display the images in layout view, similar to the graphic table elements, for each page of the map book, as shown in the image attached. I have tried the hyperlink and html popup features but they do not work in layout view. At the moment, the images you see are inserted individually but I have 100s of pictures and I would like to automate this.

One way I tried to get around this is to have the images displayed in the attribute table of the geotagged photos feature class itself and then used the graphic table element to display the corresponding fields but Im unsuccessful so far in trying to get the attribute table to display images.


texts:
I would like to get started writing a custom TagResolver, and I'm fairly unexperienced in Java.
Is there a osm2po Java API reference somewhere for the 'Default' classes (DefaultJoinInterceptor, DefaultWayTagResolver...) listing the overridable methods and their function, or some examples?

I found the TrafficSignalsNodeTagResolver.java in the plugins directory and I was wontering if there was something similar for the WayTagResolvers.

Thanks


texts: I am a beginner in GIS (for graduate school) and am using the GIS Tutorial Basic Workbook 1 for ArcGIS10.  In Chapter 1 of this text, there is a "Convert labels to annotation" section which I cannot follow.  When I right-click the item, "Convert labels to annotion" is not an option.  Please let me know how to implement this process.  Thank you for all of your help.


texts: Displacement when clip raster by shapefile

texts:
I run a script to clip raster by shapefile. Even with layers own the same crs (epsg: 4674) the output raster be displaced. I didn't find out in my script what is caused the error.

Here is my code:

"""Clip a raster image using a shapefile"""

import gdal, gdalnumeric
import shapefile
from PIL import Image, ImageDraw

def imageToArray(i):
    """
    Converts a Python Imaging Library array to a gdalnumeric image.
    """
    a=gdalnumeric.numpy.fromstring(i.tobytes(),'b')
    a.shape=i.im.size[1], i.im.size[0]
    return a

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)

def cut_raster_array(raster, shp, output):
    # Load the source data as a gdalnumeric array
    srcArray = gdalnumeric.LoadFile(raster)

    # Also load as a gdal image to get geotransform (world file) info
    srcImage = gdal.Open(raster)
    geoTrans = srcImage.GetGeoTransform()

    # Use pyshp to open the shapefile
    r = shapefile.Reader("%s" % shp)

    # Convert the layer extent to image pixel coordinates
    minX, minY, maxX, maxY = r.bbox
    ulX, ulY = world2Pixel(geoTrans, minX, maxY)
    lrX, lrY = world2Pixel(geoTrans, maxX, minY)

    # Calculate the pixel size of the new image
    pxWidth = int(lrX - ulX)
    pxHeight = int(lrY - ulY)

    # Multi-band image?
    #Check this modification in script in: http://karthur.org/2015/clipping-rasters-in-python.html
    try:
        clip = srcArray[:, ulY:lrY, ulX:lrX]

    # Nope: Must be single-band
    except IndexError:
        clip = srcArray[ulY:lrY, ulX:lrX]

    # Create a new geomatrix for the image
    geoTrans = list(geoTrans)
    geoTrans[0] = minX
    geoTrans[3] = maxY

    # Map points to pixels for drawing the county boundary
    # on a blank 8-bit, black and white, mask image.
    pixels = []
    for p in r.shape(0).points:
      pixels.append(world2Pixel(geoTrans, p[0], p[1]))
    rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
    # Create a blank image in PIL to draw the polygon.
    rasterize = ImageDraw.Draw(rasterPoly)
    rasterize.polygon(pixels, 0)
    # Convert the PIL image to a NumPy array
    mask = imageToArray(rasterPoly)

    # Clip the image using the mask
    clip = gdalnumeric.numpy.choose(mask, (clip, 0)).astype(gdalnumeric.numpy.float32)

    # Save clipping as tiff
    gdalnumeric.SaveArray(clip, "%s.tif" % output, format="GTiff", prototype=raster)
    return clip

# Raster image to clip
raster = "raster/sema/custos_it_sema_2015xmunicipio_uc_multiplicacao.tif"
# Polygon shapefile used to clip
shp = "vector/teste_corte/clip_numpyarray.shp"
# Name of clipped raster file(s)
output = "vector/teste_corte/teste_corte"

array_cropped = cut_raster_array(raster, shp, output)


texts: Is there a way to add an attribute table to a qgis2web export?

texts:
I use qgis2web for mapping building type in some area then publish it on web, so it needs a list (or table) to show on web that only a few building (not all building in that area) that I used.

Is there any way to show table (or list) of building that I used through qgis2web leaflet?

edited: I'm using qgis 2.8.3 (Wien) and qgis2web 1.17


texts: StoryMap : following a path not a location point

texts:
I am looking for a way to make a storymap with StoryMap JS (Knight Lab) or any other Opensource solution, where it can focus on a Path (a part of a waypoint), and still have the link with the text story and the next/previous step.

The solutions I have seen are  focusing on Point location, maybe I did not have a deep look.

So what I think is to adapt the Story Map JS lib get focus on path. Any one have experiences on this lib to tell me if it's hard job or not ? How to begin ?


texts: Correctly mapping CSV with x y coordinates when converting to correct CRS doesn't work

texts:
I have a CSV file that I'd like to map using lat long coordinates. The file has a field containing four coordinate points (see image below), but the first two and second two are nearly identical, so I just take the first two and separate into different columns. Then I import into ArcMap and am able to map. The points don't display properly on the map even after I convert to the CRS they were originally in.

What else can I try here?


texts: QGIS layer tree API returns only groups and not layers

texts:
I am trying to use the QGIS python API to read out the groups and layers in a QGIS (2.4) project. Here is a code snippet from my script, which is run from a shell, not from the QGIS python console!

from PyQt4.QtCore import QFileInfo
from qgis.core import *

proj=QgsProject.instance()
proj.read(QFileInfo("/var/qgis-projects/testproject.qgs"))

print "QGis version:"+QGis.QGIS_VERSION
print "Project title: "+proj.title()
print "Project tree root children:"
root=proj.layerTreeRoot()
for node in root.children():
    print "- "+node.name()+" ("+str(type(node))+")"

It prints the following when I run it:

QGis version:2.4.0-Chugiak
Project title: This is my project title
Project tree root children:
- Field Data (<class 'qgis._core.QgsLayerTreeGroup'>)
- Aerial Photos (<class 'qgis._core.QgsLayerTreeGroup'>)
- Other Maps (<class 'qgis._core.QgsLayerTreeGroup'>)
- Model Results (<class 'qgis._core.QgsLayerTreeGroup'>)

Only QgsLayerTreeGroup objects are returned but no layers (QgsLayerTreeLayer objects), even though I also have two layers in the project root. When I open the python console in QGIS while I have my project open and run:

QgsProject.instance().layerTreeRoot().children()
[<qgis._core.QgsLayerTreeGroup object at 0xa3f9170>, <qgis._core.QgsLayerTreeGroup object at 0xa3f9200>, <qgis._core.QgsLayerTreeGroup object at 0xa3f9290>, <qgis._core.QgsLayerTreeGroup object at 0xa3f9320>, <qgis._core.QgsLayerTreeLayer object at 0xa3fe320>, <qgis._core.QgsLayerTreeLayer object at 0xa3fe830>]

This does return the QgsLayerTreeLayer objects for my layers that I'm after. How come there is a difference?

Edit: I'm trying to run my python script from a shell, not from the QGIS python console.


texts: Create new polygon layers from postcodes latitude & longitude

texts:
First at all, I want to say that I am geocoding enthusiastic but have no formal background on GIS. I'll do my best to use the right wording and make myself clear what I am striving for.

Also the only tool I have available for this task is Python. This is a project I am working on in my own time but hopefully will kick off at work if I manage to produce something meaningful.

I have downloaded the Code-Point® Open from the Ordnance Survey that contains postcode units (~1.7M) for Great Britain and converted eastings & northings to Lat & Lon in WGS84.

From postcodes I've also created sector, district & area levels and a custom level that regroups areas into North, South, East & West. So my final product is a csv file that contains the levels mentioned above plus Lat & Lon. The granularity of my file is postcode level

What I want to do is to create a layer for each level and then use them later on TileMill or google map in order to create a thematic map where each patch/polygon will change color based on specific events.

In my view, what I need to do is to group these coordinates based on each level, create a collection of polygons and then create or export them as a layer (or shapefile??).

I've been looking to Shapely but I am not sure whether is the right tool for it and online examples don't make much of a sense to me  due to my lack of knowledge on this matter.

Can anyone please point me out on the right direction or provide me with a small piece of code as example?

Any advice will be highly appreciated!

[EDIT 1]

I created my layer as per the below instructions but can't see my area polygons yet.

[EDIT 2]

I've all over again 3 times and still don't get to see the polygons - I've checked my csv file and the area column has values as expected (120 areas). Here it is another screenshot with layer metadata.

And here it is a screenshot with all layers apart from original - it shows nothing.

[Edit 3]

Query option disabled.

[Edit 4]
After 4hrs processing I get to see my polygons:

In summary these are the steps:

Import csv file
Export file as shapefile
Re-import shapefile (could be new on a new project)
Apply the query logic
Apply the Convex Hull and wait.
that's pretty much


texts: Altering PostGIS query in layer data source using ArcPy with ArcGIS Pro

texts:
I've loaded several layers from a PostgreSQL database with the PostGIS extension installed into an ArcGIS Pro 2.4 project.
I'd like to make the same map a few times for the data in several different identically structured schemas.
The layer data sources are defined using PostgreSQL queries.
for example I need to change SELECT * FROM california.scratch.d4_level_1_needs_pt
into SELECT * FROM california.scratch.d5_level_1_needs_pt
or a slightly more complicated query:
SELECT
    geom,linkid,route,finalscore,urgency,mobilityscore,
FROM
    california.scratch.d4_level_1_needs_ln

It looks like I can update the layer source using:
import arcpy, pprint
p = arcpy.mp.ArcGISProject('current')
m = p.listMaps()[0]
l = m.listLayers()[0]
pprint.pprint(l.connectionProperties)

and arcpy.mp.ArcGISProject().updateConnectionProperties()
but when I print the connectionProperties of one of the layers I don't see the query. it gives me:
{'connection_info': {'authentication_mode': 'DBMS',
                     'database': 'california',
                     'db_connection_properties': '***.***.**.**',
                     'dbclient': 'postgresql',
                     'instance': 'sde:postgresql:***.***.**.**',
                     'password': '<*********>',
                     'server': '***.***.**.**',
                     'user': 'user'},
 'dataset': 'california.scratch.%d4_partners_surveys_3',
 'workspace_factory': 'SDE'}

I could update the value for key dataset but that doesn't help me if I want to update the columns I'm pulling in via the query.
Additionally, is there any way to associate the connections I print in the python window with the layers in the list by Drawing Order?
I don't see any unique identifying info.


texts: Delete duplicate spatial point/points in mssql from a perticular location

texts:
I am working with a Spatial data which has a street intersection points....
I get the intersection points from spatial query on street/roads datain MMSQL . but in my case some streets/roads are multi lane.. which gives me duplicate intersection points for same streets/roads for each lane about 3-5 meters away from each other..
now i have the point data with Duplicate intersection points with a geometry differing in about 3-5 meters.. Now i want to delete the duplicate intersection points and keep only one intersection point.. i tried taking buffer of a point n den deleting the points which intersects with this buffer.. the the query is very slow and takes ages to complete.. i tried using spatial indexes but it hardly helps...

is there a algorithm or any method in QGIS ARCGIS OR POSTGIS Database which i can use to delete theses duplicate points.. or any other method?????


texts: Cannot export using "Save Selected Feature As..." in QGIS because the name is more than 10digits. How to export without manually changing the column

texts: I am trying to export the shapefile that I have to a smaller shapefile using "Save Selected Feature As...". But, the saving process is error due to the length name of the column. How can I export it and get the shapefile without manually changing the column attribute name manually?


texts: How to find and replace attribute values based on a list

texts:
I have a Shapefile and a list in .csv which contains two columns with names:

"Old_name" which contains the old names that already exist in the attribute table of the previously mentioned Shapefile; and
"New_name" which contains the updated name.

How can I do this either on ArcGIS, QGIS or with python?

Based on this question I need to use the UpdateCursor function, but how can I use the list in this function?


texts: How to load a style from POSTGIS in PyQt

texts:
The thing is that I have four styles I have created in QGIS through the GUI and saved into a postgresql database.

I can see them in the table named layer_styles

I have seen similar questions asking to load from qml file, or even from a database, like that but it looks as if I need to get a styled layer before for get the style from it.

It's like I say or can I load the styles from database without a styled layer?

Anyway I don't know how must I do to get the styles and use it for another layers


texts: Count number of attributes for unique field values in ArcGIS for Desktop?

texts:
Running ArcGIS 10.1 for Desktop, I have two feature classes in a geodatabase.  One is a point FC for polling stations, the other is a polygon FC for states.

The polling station FC contains over 1000 points and has a states field indicating which state the polling station is located in.  The states FC has an attribute for each state and has a blank field for the number of polling stations in each state.

I need a python script that will get the count of how many polling stations are located in each state from the polling station FC, and then write those numbers to the column in the states FC - obviously to their corresponding state.


texts: I want to style a point layer from png file. Can i do this in quantum gis?


texts: Finding nearest points to POI with big data at Yelp?

texts:
Businesses often need to find nearest points to Points of Interest (POI).

For example, use Yelp to find restaurants within 10 mile distance. As a normal GIS practice, we solve this problem by using the GIS functions. There are lots of Q&As here on this topic and examples are:

Finding points near points using PostgreSQL?
Finding nearest point from POI in PostGIS?

I wonder how about in the real business world, when data cannot fit in memory and high computational efficiency is a must.

How does Yelp solve these kind of spatial problem in their work?

Perhaps they use PostGIS, QGIS functions, etc to find nearest points near POI, or maybe they develop their own tools.


texts: Is Shooting Star broken in PGRouting 1.05?

texts:
I'm using pgrouting DLL 1.05 with postgresql 9.1.
All required function installed.

Astar and Dijkstra work fine with sample data but "SHOOTING STAR" is not considering "reverse_cost".

Also if i try long route then it gives server disconnected error.

"Do you want to attempt to reconnect to the database?"

And then I need to restart the postgresql service.


texts: how to search for a Point of MultiLineString geometry in a json file using OpenLayers

texts:
I have a json file with large geometries(MultiLineString), here I got a Point among the json layer(baseLayer). I would like to search for a Point's LineString in the layer, from which MultiLineString that Point been selected. And how to read a json file in OpenLayers/JavaScript

I can project the json file in OpenLayers, but when I click the map I am getting the Point, from that Point I would like to search for its base geometry(LineString/MultiLineString)

Thank you


texts: Vector File Merging fails as string length too long

texts:
I am running QGIS 3.18.  I am finding my feet with QGIS.  I am building up maps of my forebears' farms. From death duty and probate documents I have collected the lot numbers of the paddocks. I have then downloaded the corresponding shapefiles from our state lands department.
With properties with acreages in their thousands, I have lots of shapefiles to manage.  Using Vector>Data Management Tools> Merge Vector Layer, I created a file (.qpkg) which represents the freehold acreage for a particular year.  I then tried doing this for the lots representing the leasehold holdings and have run into a problem - field PARCEL_SPI is defined as a string of 10 characters, whereas about 2/3 of the lots I want to merge for the leasehold have a PARCEL_SPI string field containing 11 chars. (whether the lot is leasehold or freehold is irrelevant - I was lucky with the freehold that all of the PARCEL_SPI values were either 9 or 10 characters long.)  The data in PARCEL_SPI is of the format of Lot Number\PPxxxx. This is fine for the case of 34A\PP2846, however gives errors in the case of 118A\PP3618 because it is 11 characters long and the vector file merging process has a field size of 10.
Presumably the easiest way is to define the Qstring for PARCEL_SPI as 11 characters, how do I do this?


texts: Agisoft Metashape Orthomosaics. DEM from mesh or dense point cloud

texts: Does anyone know the differences between Orthos built from DEMs that were constructed from dense point clouds rather than mesh? My purpose is to preserve original data and minimize interpolation while having a continuous ortho.


texts: I'm getting started with a qgis map showing the quality of roads in my town. I have various data which create an index of 1-5. How do I link this index data to road data in a map, showing the bad roads (5) in red, and good roads (1) in green? What is the correct term for this process? Any pointers on this would be greatly appreciated. Regards, Case.


texts:
Is the line color and thickness of the ArcGIS measure tool able to be customized?

My searches have come up empty.
\
I'm working against a dark background and the black line (from the measuring tool) does not provide enough contrast to distinguish between the measure tool and the background.


texts: Editing / grading a DEM raster? (Ex: detention pond and berm grading)

texts: I am looking for a way to adjust/edit the values in a DEM file, for example, to grade a detention pond and berm into existing LiDAR data. I have tried messing around with point shapefiles and the IDW tool, but I am wondering if there is a better way. I can provide any additional details as needed.


texts: Distortion of an exported image in Google Earth Engine

texts:
I am currently analysing the MODIS datasets (e.g "MODIS/006/MOD16A2") available in Earth Engine, computing mean values for each year. While the original projection of the MODIS datasets is SR-ORG:6974, I understand that when functions like mean() are applied to the image collection, the projection changes to EPSG:4326. This transformation does not concern me greatly as I would then apply a reduceRegions() function over the image using a Feature Collection (grid layer) that was based off the EPGS:4326 projection; aggregating the values within the 1km grid squares.
After converting the Vector (from the reduceRegions) back into an Image, I exported the image into my assets. However, when analysing the exported image, I realised that the exported image does not match the computed image in the Map viewer, resulting in incorrect data input for that grid square.
Hence, I would like to figure out if the issue arises based on the projections of the computations or could it be an issue with the export of data.

From the image, that there is an additional layer of white grid squares that were created. This has led to those grid squares holding different values than the original one.
Line of code used:
//------------------ Extracting THA Boundary -----------------------------

//  Call ee.FeatureCollection (i.e. vector data)
var ssea = ee.FeatureCollection('users/jjohanness1992/ssea/ssea-adm0');

//  Filter ee.FeatureCollection only to Southeast Asia
var countrylist = ee.List(['THA']);
var THA = ssea.filter(ee.Filter.inList('GID_0',countrylist)); // Filters to Thailand

Map.centerObject(THA,6);

//-----------------------------------------------------------------------

var evapotranspirationVis1 = {
  min: 0.0,
  max: 300.0,
  palette: [
    'ffffff', 'fcd163', '99b718', '66a000', '3e8601', '207401', '056201',
    '004c00', '011301'
  ],
};

var evapotranspirationVis2 = {
  min: 0.0,
  max: 1600.0,
  palette: [
    'ffffff', 'fcd163', '99b718', '66a000', '3e8601', '207401', '056201',
    '004c00', '011301'
  ],
};

var evapotranspirationVis3 = {
  min: 0.0,
  max: 5500.0,
  palette: [
    'ffffff', 'fcd163', '99b718', '66a000', '3e8601', '207401', '056201',
    '004c00', '011301'
  ],
};

var ET = ee.ImageCollection("MODIS/006/MOD16A2")
  .filterDate('2018-01-01','2019-01-01')
  .filterBounds(THA)
  .select('ET')
  .map(function(image) { return image.clip(THA)});

print('ET',ET);
print('Projection of ET:', ET.first().projection());
print('Scale of ET:', ET.first().projection().nominalScale());

Map.addLayer(ET,evapotranspirationVis1,'ET');

var totalET = ET.sum().multiply(0.1); // adjusted for scaling factor
print('Projection of Total ET:',totalET.projection());
print('Scale of Total ET:',totalET.projection().nominalScale());

Map.addLayer(totalET,evapotranspirationVis2,'annual total ET');

// summing up into 1km size
var tha_grid = table;

var stats = totalET.reduceRegions({
  reducer:ee.Reducer.sum(),
  collection:tha_grid,
  scale:463.3127165275
});

var totalET1km = stats.reduceToImage({
  properties: ['sum'],
  reducer: ee.Reducer.first()
});

Map.addLayer(totalET1km,evapotranspirationVis3,'ET_1km');

Export.image.toAsset({
  image:totalET1km,
  description:'THA_TotalAnnualET_2018',
  assetId:'Thailand/ET_Sum/THA_TotalAnnualET_2018',
  scale:926.6254330555,
  region:geometry
//Here is a geometry that I made myself which is a box over the Thailand region. I've also tried using
//THA.geometry().bounds() but I also get the same results
});


texts: Adjusting extent of raster using ArcGIS Desktop?

texts:
The coloured file contains landuse information and it needs to have the same number of columns and rows as the DEM in the background for further processing.
I can't get it to work - "clip" and "extract by mask" haven't generated the results I need.

It's supposed to keep its values within the coloured areas - but outside of that to the borders of the DEM (the whole black and white area) I need it to have NoData values.


texts: Display Features based on Zoom Level using OpenLayers

texts:
I want to display roads(linestring) on the map using OpenLayers.currently, I manage to assign styling based on road type but when I run the server I can see all the data displayed. I want to change the way it displays feature.
for e.g.

CATA roads(orange color) should display above zoom level 12,
CATB roads(yellow color) should display above zoom level 14 but at the same time, it should also display CATA roads.
Other roads should display above 16 but it should also show CATA and CATB roads along with it.

In the current code I am creating a function to find the road category and assigning style to each category.
I tried to incorporate zoom level and minZoom in my code but I am not able to make it work. It doesn't change anything the way it displays.

code:
var styleRoadFunction = function(feature,resolution) {
      var zoom = map.getView().getResolution();
      var road_categ = feature.get("road_categ");
      //var zoom = map.getView().getZoomForResolution();
      console.log(feature);

        var color;
        var style;
        var text;
        var minZoom;

        //if (zoom >=12){
        if (road_categ == 'CATA' && zoom>=11.9999 ){

          style = [
            //for outer color
          new ol.style.Style({
            stroke: new ol.style.Stroke({
              color: '#F6B25C',
              width: 10,
              lineCap: 'butt',
              lineJoin:'round',
              zIndex:1,
              //minZoom: 12,
            })
          }),
          //for inner color
          new ol.style.Style({
            stroke: new ol.style.Stroke({
              color: '#F6B35C',
              width: 7,
              lineCap: 'butt',
              lineJoin:'round',
              zIndex:1,
              //minZoom: 11.9999,
            })
          }),
        ]
       }
        else if (road_categ == 'CATB' && zoom>=13.9999){

          style = [
              //for outer color
            new ol.style.Style({
              stroke: new ol.style.Stroke({
                color: '#DED5B5',
                width: 9,
                lineCap: 'butt',
                lineJoin:'round',
                zIndex:1,
                //minZoom: 12.9999,
              })
            }),
            //for inner color
            new ol.style.Style({
              stroke: new ol.style.Stroke({
                color: '#F9E7A1',
                width: 7,
                lineCap: 'butt',
                lineJoin:'round',
                zIndex:1,
                //minZoom: 12.9999,
              })
            }),
          ]
        }
        else
          style = [
          new ol.style.Style({
            stroke: new ol.style.Stroke({
              color: '#c4c2be',
              width: 7,
              lineCap: 'butt',
              lineJoin:'round',
              zIndex:0,
            })
          }),
          new ol.style.Style({
            stroke: new ol.style.Stroke({
              color: '#ffffff',
              width: 5,
              lineCap: 'butt',
              lineJoin:'round',
              zIndex:0,

            })
          }),

          return style;
         };
  var vectorLayer6 = new ol.layer.VectorTile({
            source: new ol.source.VectorTile({

            format: new ol.format.MVT(),
            url: 'http://localhost:8081/tiles/roads/{z}/{x}/{y}.mvt'
        }),

        style:styleRoadFunction
       })


texts: One-liner to print first feature from GeoJSON

texts:
To make the findings from this GDAL issue more discoverable, I'm posting here as well.

ogr2ogr has a -limit option. I couldn't find this option for ogrinfo, while ogrinfo has e.g. -where. I would like to use the same simplicity like with the -where flag, e.g. neither specify select * nor from layer.
I use it to look at actual values of a layer. -so only gives the column names and column types, and simply ogrinfo spams the terminal if you are doing it on a file with 2.5 mio features.
I'm working inside a JupyterLab and there is no way to start any GUI applications. There are only a Terminal, iPython console and Jupyter notebooks. Besides, Ctrl+C to stop the spamming might not be carried out for quite some time.
Simply doing head file.geojson does not work if the file has all its contents on one line.
How can I do something like ogrinfo -limit 1 (ogrinfo has no -limit option) for a GeoJSON file?


texts: Adding shapefile or feature class as layer in ArcMap using ArcPy

texts:
I am trying to automate various tasks in ArcGIS Desktop (using ArcMap generally) with Python, and I keep needing a way to add a shapefile to the current map. (And then do stuff to it, but that's another story).

The best I can do so far is to add a layer file to the current map, using the following ("addLayer" is a layer file object):

def AddLayerFromLayerFile(addLayer):
 import arcpy
 mxd = arcpy.mapping.MapDocument("CURRENT")
 df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
 arcpy.mapping.AddLayer(df, addLayer, "AUTO_ARRANGE")
 arcpy.RefreshActiveView()
 arcpy.RefreshTOC()
 del mxd, df, addLayer

However, my raw data is always going be shapefiles, so I need to be able to open them. (Equivalently: convert a shapefile to a layer file without opening it, but I'd prefer not to do that).


texts: ArcPy SearchCursor Loop Only Executes Twice

texts:
I'm trying to loop through a set of lines and select points within .25 miles of them. I just need the count, so I'm trying to put the results into a list I can print out. The relevant code is below.

When I execute the code, it will loop through twice and then print the list called data. There are actually 35 rows it should loop through. If I comment out everything but printing out the name of the route (first line of the loop), it will successfully execute as I intended.

I found a similar question here: Select Layer By Attribute using loop. However, this question simply ended with the poster figuring out that their cursor was tied up. Not sure exactly what that means or how to resolve it.

Why would the loop execute less times than there are rows, but not produce an error?

fc = r'PATH TO FILE'
inlyr = r'PATH TO FILE'
data = {}
field = 'ROUTE'

arcpy.MakeFeatureLayer_management(fc, 'route')
route = 'route'

arcpy.MakeFeatureLayer_management(inlyr, 'addrs')
addrs = 'addrs'

cursor = arcpy.SearchCursor(fc)

for row in cursor:
    print "ROUTE = \'{}\'".format(row.getValue(field))
    arcpy.SelectLayerByAttribute_management(route, "NEW_SELECTION", "ROUTE = \'{}\'".format(row.getValue(field)))
    arcpy.SelectLayerByLocation_management(addrs, 'WITHIN_A_DISTANCE', route, '.25 Miles')
    num = arcpy.GetCount_management(addrs)
    data.update({row.getValue(field):num})
    del num
del cursor
print data


texts: How does QGIS know there is a plugin update?

texts:
I've created a plugin and hosting the installation zip and xml file on a local network repository (for company use). Everything works fine but every time a user opens QGIS, it shows that the plugin has an update when nothing has been updated. The issue looks like this:

Is there something in the repository xml file that flags QGIS that the plugin has an update available? How can I control this so it only flags when there is actually an update?


texts: Large loads geojson from mysql to leaflet?

texts:
I have a mysql database with +25K rows including geometry and other attributes.
Using a mysql2geojson php script I then added it to a leaflet map.

It is super slow now, what can I do to boost the performance?

I'm using www.one.com as host meaning I dont have access to anything else than php/mysql. I can't get postgis/geoserver, etc.

I also intend to design a search filter limiting which points to be shown on the map. For this purpose, is it best to do new calls to the database, or should I be working against a "local dump" for the searching?

To summarize I wonder how to best work with big point data sets with mysql/php/leaflet including searching. Any tips are welcome!


texts: NDVI Analysis of a small area from the year 2000- 2011

texts:
I want to perform an NDVI analysis in a small area. I want to see the change around a river body for the year 2002, 2005, 2008, 2010 and 2011. I tried using Landsat-7 data but due to the scanline error, it's not a good choice.

I want to work on free data for now.

Which satellite data should I use?


texts: I am trying to run TauDEM algorithims with zero luck. I am running QGIS 2.4, have installed TauDEM 5.0.6, MPICH2 and still am getting nowhere. I either get a problem loading output layers or directory name is invalid. Please help.


texts: how use raster data directly with openlayers

texts: I am currently using qgis2web to set up a WebGIS project for which I would like to avoid a WebGIS server. I cannot find a way to convert my raster data in a way that it keeps its information (1 band 32 bit float). If I use the qgis2Web export function, it gets converted into a .png without the original values.  I tried to search a file format that is supported by Openlayers or Leaflet that would retain the original values but could only find the option to use WMS.


texts: How to get an older version of OpenStreetMap data?

texts:
For some reason in the map of the Dominican Republic's Capital there are a bunch of streets which doesn't have any names. I had an OSM file which I downloaded ealier this year which did have those streets name.

I really don't know WHY the names has just dissapeared, but I would like to know if there is a way that I can get them back.

Maybe I can download an older version and then re-upload it so the current version get the streets' names back.

So can you tell me how to get an older version of the OpenStreetMap or can you give me any other alternatives to get the streets' names back?

Thanks in advance.


texts: Determining ArcGIS license server license is coming from?

texts:
I have a esri license server setup with a failover (as described in the ArcGIS documentation).

All the ArcGIS clients have the environment variable ARCGIS_LICENSE_FILE set to 27001@mainlicenseserver;27001@failoverlicenseserver.

Is there a way a can see from which license server a ArcGIS client (for instance ArcMap) got the license?


texts: How do I change geotiff "no data" color to white/transparent rather than black?

texts:
I am preparing data for inclusion into a geonetwork node, where the maps are displayed via geoserver.

The problem is that cells that have"no data" are displaying as black in the geotiff and I am would like the color to be set to white or transparent.

I am using ArcView and do not have access to either Spatial Analyst or 3D analyst.

I would prefer not to have to assign a new sld file and am presently using "Use colormap" when exporting from ArcGIS Desktop.


texts: Filtering layers that show for input feature layer parameter for custom tool?

texts:
When creating a custom tool using a Python script, is there a way to only show layers that, for example, contain a certain word?

For the drop down of the input polygon feature, I only want layers to show that have the word boundary in them.

Is this possible, or does it just automatically show all layers in the document?

Here's a screenshot of what I'm talking about:

It's a bad example because it only shows one layer, but image lots, and lots of layers.


texts: Connect to Enterprise GDB in PostgreSQL using QGIS : "No Geometry" Error

texts:
I have an esri Enterprise Geodatabase (9.3) stored in MSSql 2008 with storage type as sde binary.

I upgraded it to EGDB 10.4.1 (using multiple up-gradation steps) in MSSql 2012.

Then I copy/paste my data to PostgreSQL(9.4). I have the storage type in PostgreSQL as high precision ST_Geometry spatial type (SRID 300001).

I tried to connect to PostgreSQL geodatabase using QGIS but all my data are showing as no geometry and hence I cannot view it in QGIS.

However I can connect to it using ArcMap.

How do I solve this?


texts: How to lock the legend in print composer?

texts:
I locked in the print composition (Druckzusammenstellung) the legend (attachment1).
Then I saved the print composition and opened it (attachment2). But the locked legend is not locked? The changes are gone: polygon b is still part of the legend.
What have I to do to lock the legend?

Thank's for any advice.
phil

attachment1:

attachment2:


texts: How to show Ocean Contours on Mapbox Style?

texts:
My users would like to be able to see the ocean contour lines as they zoom in. But mapbox outdoors style only seems to provide contours for land.

I'm not a GIS pro but a web dev, so please keep that in mind in any answers. I do have experience in publishing my own tilesets to mapbox with points, shapes, and labels. But a tileset with all of the contours for the entire water area of the earth sounds impractical and very large for me to host on my little portion of mapbox. But I don't see that they have such a layers available.


texts:
We are currently running ArcGIS 10 in our company through Citrix as well as installs on user workstation.  We always seem to have niggling problems with Citrix hence local installs.

For example:

Poor display of smoothed grids in Citrix showing very blurry but then perfect display when exported to a PNG
Mapped drives not being able to connect but once the session is closed down and re-opened and it's fine again
ArcScene has terrible performance when it has a couple of stacked surfaces to render
ArcGIS Extension licenses not always being released when switched off in extension dialog (user has to shut down their session)

There are lots more issues that can be found on the forums of course and I have read the strategy documents too but I'm looking for more detail on how to actually deploy it without these niggling problems!

But obviously as a concept it is great as it saves so much time for IT and managing the application centrally is the way forward, etc. But the users are suffering at the moment.

So my question is what are peoples opinions with deploying ArcGIS through Citrix?

Have people had "strategists"/Citrix experts to come in and show us how to set it up specifically for ArcGIS to get the best performance and display?

Basically the users should not see a difference between ArcGIS on their machines and ArcGIS on Citrix. Otherwise it won't get used.


texts: Creating new layer from selection in ArcPy?

texts:
arcpy.SelectLayerByAttribute_management() helps in selection but I am unable to create a new layer from the selection arcpy.MakeFeatureLayer_management() for generating new layer file but it gives an error:

arcpy.MakeFeatureLayer_management("C:/Users/Gurminder/Documents/study/IWMI/TestingAutomation/SnowLine/SnowPersonal.mdb/Snow_lyr")
Runtime error  Traceback (most recent call last):   File "<string>", line 1, in <module>   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 5774, in MakeFeatureLayer     raise e ExecuteError: ERROR 000732: Input Features: Dataset C:/Users/Gurminder/Documents/study/IWMI/TestingAutomation/SnowLine/SnowPersonal.mdb/Snow_lyr does not exist or is not supported ERROR 000735: Output Layer: Value is required

How do i create a new layer so that i can copy the selection in the newly created layer, I have used Personal geodatabase

Is there any alternate way to do the same?


texts: rscript in QGIS processing toolbox and Rstudio, results mismatch

texts:
I am using my own R scripts in QGIS. Before adapting them to work in QGIS, I write them in RStudio to  make sure they work properly. I then use different scripts in a batch model. In RStudio I have loaded the same files from the same directories that I have also loaded in QGIS. With my code I load shapefiles and rasterfiles, make some operations, save some output shapefiles and rasterfiles and then automatically produce a TIFF image and save it. The filename of TIFF includes part of the original filename of one of the shapefiles loaded into R. When I run my code in Rstudio, everything works as expected. When adapting the code to QGIS, the TIFF does not have part of the original filename of the desired shapefile, instead it is replaced by a "NA". The whole script produces otherwise the expected results. For the love of money I cannot figure what goes wrong there. My code is:
#load my files
shp<-readOGR("myshapefile.shp")
rast<-raster("myrasterfile.TIFF")

#run my calculations
#save the desired output shapefiles/rasterfiles

#retrieve the original filename of the shapefile loaded as "shp"

    ab<-as.data.frame(shp)
    abc<-rrapply(as.list(.GlobalEnv), classes = "data.frame", how = "flatten", dfaslist = FALSE) #from package rrapply which as far as I understand does not work with S4 class objects, that's why of previous command

    outfile <-paste0("NewName_", names(abc[1]), ".tiff") #I have also tried with other numbers for names(abc[]) just to make sure I wasn't searching in the wrong place, but nothing changes

#set the TIFF parameters
tiff(filename=outfile,
res = 800,
compression = "lzw",
width = 5,
height = 5,
units = "in"
)

#make the plot

The saved TIFF should be named "NewName_myshapefile.tiff" and it is instead named "NewName_NA.tiff" and this is because the variable abc is empty, in fact, if I run
names(abc)

I get an NA as result. In R I do not have this issue and everything goes as expected. I am 100% positive that my model is not the source of error, if I run my script alone I get the same result.
How do I solve this?
EDIT
based on my comment to @Spacedman, here is my new code to generate the desired names for the output TIFF files:
packages <- c("rgdal", "adehabitatHR",  "raster", "rasterVis", "sp", "sf", "lattice", "viridis", "spData")

package.check <- lapply(
packages,
FUN = function(x) {
if (!require(x, character.only = TRUE)) {
install.packages(x, dependencies = TRUE)
library(x, character.only = TRUE)
}
}
)

#make lists of raster and vector layers from which Iwill be matching the name to generate the output TIFF filename
rlist = list.files(getwd(), pattern = "tif$", full.names = FALSE)
shplist = list.files(getwd(), pattern = "shp$", full.names = FALSE)

#run my calculations
#save the desired output shapefiles/rasterfiles

#retrieve the original filename
numb<-grep("mypatternforthenamegoeshere", shplist)
outfile <-paste0("Newmane_", as.character(shplist[numb]), ".tiff")

tiff(filename=outfile,
res = 800,
compression = "lzw",
width = 5,
height = 5,
units = "in"
)

#make my plot

When I adapt the script to work in QGIS, amongst my input files I add a string character field where I type the pattern I want to match in the list, job done. In QGis syntax this is done by declaring at the beginning of the script:
##Make_name_based_on=string

Therefore, in my edited code above, the line:
numb<-grep("mypatternforthenamegoeshere", shplist)

will become:
numb<-grep(Make_name_based_on, shplist)

This gives me what I needed.


texts: How to insert new data points into CartoDB map?

texts:
I'm using NASA api data to populate my database in CartoDB's website to load places where meteorites have fallen.
In JavaScript I'm currently creating a specific collection of meteorites I want to appear on the map, and these will change based on how users interact with the map. The more they click, new ones should appear.

The issue I'm facing currently is that I cannot find out how to call for this specific collection of meteorites to appear on the map. Since it changes based on several factors, I'm thinking that making SQL queries won't work and instead I would want to create a new viz.json file.

I'm confused on where does the data go into cartodb.createviz(. Looking at examples such as visjson spec there's no hints of where exactly the data goes.


texts: Gdal Read Geodatabase Features Don't Show Object ID

texts:
I am using Gdal Node.js interface to read features of a layer from a Esri Geodatabase. My javascript code looks like this:

var dataset = gdal.open(inputPath);
var layer = dataset.layers.get(layerName);
layer.features.forEach(function (feature, featureId) {
    var geom;
    var recordId = null;
    var routeId = null;
    try {
        var geomWKB = feature.getGeometry().toWKB();
        var geomWKBString = geomWKB.toString("hex");
        geom = geomWKBString;
        if (recordColumn) {
            let names = feature.fields.getNames();
            recordId = feature.fields.get(recordColumn);
        }
       if (routeIdColumn) {
           routeId = feature.fields.get(routeIdColumn);
       }
    }
    catch (err) {
    }
    ......
});

This layer has 12 fields. However, gdal only gets 11 fields missing OBJECTID when extract fields of a feature. Therefore recordId = feature.fields.get('OBJECTID') does not work.

I inspected the layer and noticed OBJECTID was assigned to FIDColumn, which means gdal can access to OBJECTID field but doesn't add it to features' schema.

What could cause the issue?


texts: Apportioning Points at a street intersection (arcgis or mapinfo)

texts:
I can't seem to find any evidence of this question, but I'm surprised by that, so pardon if this has been asked.

If one is mapping crime location reports, by XY coordinates to streets, to create "hot street" segments, how does ArcGIS handle intersections? Which street (the East/west or north/south) and/or which side of the street, gets counts of events at intersections?

Is there a default logic to this? Looking at my map some segments get both intersections, some get one, and they all go to east/west streets, not the perpendicular intersecting streets.

For the left/right (east/west) issue, where some segments get both intersections and some get one, my guess was it's based on a true north/south/east/west grid, and whether the intersection is to the east or west, but I can't figure out why the perpendicular streets don't seem to get intersection incidents.

Maybe a work around is having those counts be specific to the intersection and not to the line segment? But, I'm not sure how to do that.


texts: arcgis geodatabase replica between different versions

texts:
We have 2 geodatabase corresponding to 2 systems. There is a replica between these 2 geodatabase all use version 10. Recently there is a version upgrade in the parent geodatabase.
It will be upgraded to version 10.2.2 (the child veision will remain version 10). I want to confirm can I create a replica between version 10 and version 10.2.2. Is there any potential danger?

We have a scheduler to synchronize changes of the 2 geodatabase use the replica now. with the upgrade of the parent geodatabase, will the replica still work?


texts: creating multi band raster from one GeoTiff

texts:
I'm using AgiSoft photoscan to create a photo mosaic from aerial images. The images are being captured using a camera modified for NIR and RGB data capture.

I want to be able to perform spectral mixture analysis on the image using a package like ENVI but am unsure where to start or if this is even possible.

I assume that I would need to first create a mutli band raster from the single exported ortho GeoTIFF from Agisoft Photoscan. But even this I'm not sure about. Can anyone point me in the right direction please?


texts: QGIS Orfeo Toolbox: Smoothing tools generate unwanted black border on a hillshade raster

texts:
I'm using the Smoothing (gaussian) and Smoothing (mean) tools in the QGIS Orfeo toolbox to blur a hillshade raster I've generated. However, they both come out with a thick black border at the edge of the DEM (essentially at sea level as well as at the edge of the DEM where I clipped it).

What I'd like to know is:

A) Can I prevent this somehow?
B) If not, is there a tool I could use to easily remove it?

I've also tried Smoothing (anidif) but it causes a lot of unwanted artifacts within the data - it causes the raster to appear to be divided into slightly discontinuous horizontal strips.

Screenshots: Raw input, and Smoothing (gaussian) output respectively. The black border with smoothing (mean) is a lot thinner, but the results of that are also inferior to the gaussian smoothing.


texts: I'd like to create flat buffer in a polygon shapefile. I try to use buffer tool in ArcGIS 9.3, but I couldn't select the "flat" option in the endtype combo-box because the input feature is polygon. How can I select a flat buffer in buffer tool?


texts: Why does QGIS QuickMapServices' export of hi-res aerial maps return blank images?

texts:
Using the QuickMapServices plugin in QGIS 3.12.0 to load Esri World Imagery (Clarity) Beta, I noticed that every time I try to export large scale maps (e.g. 1:1,000) as high-resolution images  (at least 300dpi) the result is a blank image.

At 300 dpi, the largest scale export that works fine is 1:3,000. So there are no issues with aerial images at this scale and at scales lower than that (1:5,000, 1:10,000...).

I have tried exporting various file types, but the result is the same.

My colleague, who uses the same QGIS version and the same computer configuration (with 8GB RAM), has not experienced this issue.

I have re-installed QGIS, now using the latest version (3.12.1), but it still exports aerials as blank images.

Edit (further info):

To save an aerial image in QGIS I simply go to Project -> Import/Export -> Export Map to Image
It seems that exporting other aerial images (Bing) works fine.


texts: ArcMap Modelbuilder Iterator (tables) - Not iterating DBF files [SOLVED]

texts:
I have a model created in ArcMAP 10.8 modelbuilder that has an iterator (tables) (standard setup). The model is setup in such a way that it runs various analysis on multiple DBF files in a folder. The model has always worked but now it ignores all the DBF files in the folder as if its empty.
When the model has finished I get this message "NOT_RECURSIVE"
Please advice


texts: Eliminating a column from a copied feature class via a fieldmap object in ArcGIS Python

texts:
The use case is exporting a feature class via a Python script (arcgisscripting v9.3x.  Only a few of the columns are needed in the result feature class.

I commonly use the fieldmappings/fieldmap objects to rename columns. The docs state that you can also use them to delete a field as the feature class is being exported.  Can anyone provide an example of this?


texts: Error while opening shapefile in QGIS on Surface Pro

texts:
I downloaded QGIS on Surface Pro, but when I try open a .shp-file the following
error comes up:

Layer is not valid: The layer C:/Users/Izzyg/Documents/lab
1/__MACOSX/GEOG2480_F20_Lab01Data/countries/ne_10m_admin_0_countries.shp is not a valid layer and can not be added to the map.

Reason: I have tried on another computer and I am able to open the file just fine on there. What can be a reason?


texts: Interacting between Openlayers and Python

texts:
I have a openlayers website like this one that gives the user the ability to draw polygons and then creates a GeoJSON file. I want to process this GeoJSON with a python script. The python script returns in the end a string that should be returned to the website and displayed as a info box.

What kind of frameworks/apps/steps should I look into to accomplish this project?


texts:
Is there a way to get the type of all the layers in an OpenLayers maps object? What I am trying to do is iterate over all the layers in my OpenLayers map object and populate a drop-down list with the names of only the vector layers. All my vector layers come from a remote GeoServer and are WFS.

Here is the code I currently have that gets the names of the layers and populates a drop-down with the names. This code works fine;

function listLayers() {
var mLayers = map.layers;
//set a variable to the layers array in the map object
document.getElementById("drpLayers").options.length = 0;
// clear any existing options from the drop-down list

for (var a = 0; a < mLayers.length; a++) {
    //loop through the map layers array
    var opt = document.createElement("option");
    //create an option object in the drop-down
    opt.text = mLayers[a].name;
    //add the name of the layer for this drop-down list option to the text attribute of the option object
    opt.value = mLayers[a].name;
    //add the name of the layer for this drop-down list option to the value attribute of the option object
    document.getElementById("drpLayers").options.add(opt);
    //add the option object to the drop-down list
};

I thought I could just modify this code to something like this:

opt.text = mLayers[a].type

or perhaps this:

opt.text = mLayers[a].class

but neither way works. When I look at the OpenLayers documentation I cannot find a layer property that describes the layer type. Any suggestions would be appreciated!


texts: How can I implement Koppen's climate classification in ArcGIS?

texts:
I am about to start a project to make a map showing the climate zones in US. I am given a dataset of the climate, but my next step is to make a separate layer of each of the climate groups defined by Koppen.
For example, Group B: dry climates is defined by little precipitation, and according to Wikipedia, it can be calculated by completing the following steps.
Multiply the average annual temperature in Celsius by 20, then add
(a) 280 if 70% or more of the total precipitation is in the spring and summer months (April–September in the Northern Hemisphere, or October–March in the Southern), or
(b) 140 if 30%–70% of the total precipitation is received during the spring and summer, or
(c) 0 if less than 30% of the total precipitation is received during the spring and summer.
Assume that I already have all the data I need for this calculation such as the amount of precipitation, what tool in ArcGIS would I need to complete this calculation?


texts:
I have designed a polygon in a web page and the coordinates are stored in the database of the type polygon. However when I check the database, I find GEOMETRY - 120B under polygon field.

What does it mean and how can I know the coordinates? In the end, I want to know whether a point is inside or outside the polygon?


texts:
Working in QGIS 3.4.1 Madeira

I have a layer with many data points, for visual purposes I wish to remove some of the overlaying points.

Is there a simple way to perform some sort of clustering/filtering of point based on their relative distance to each other?

I found the Plugin calls "ClusterPoint" but cannot find where it is or how to use it.


texts: Exporting a GIF of Sentinel-5P image collection with country outlines

texts:
I've been stuck on this for a couple of days, and I feel like I need a second set of eyes for this.
I have made a script that calls the S5P image collection for NOx emissions over period of time, and export a GIF for a time lapse.
My next step would be to add in country outlines, to give the viewer some spatial reference of what they are looking at.

Is there any documentation of how to do this for a GIF export?

Current Code:

*/ I used the geometry import function to draw out my area of interest.

var startDate = '2019-04-29';
var endDate = '2020-04-29';
var weekDifference = ee.Date(startDate).advance(1, 'week').millis().subtract(ee.Date(startDate).millis());
var listMap = ee.List.sequence(ee.Date(startDate).millis(), ee.Date(endDate).millis(), weekDifference);

function getWeekly(date) {
  var date = ee.Date(date);
  var sentinel1 = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2')
                      .select('tropospheric_NO2_column_number_density')
                      .filterDate(date, date.advance(1, 'week'))
                      .mean();
  return sentinel1;
}

var image_prior = ee.ImageCollection.fromImages(listMap.map(function(dateMillis){
  var date = ee.Date(dateMillis);
  return getWeekly(date);
}));
var Nox = image_prior.combine(image_prior).combine(image_prior)
.map(function(image) {
  return image.divide(0.0003).multiply(255).uint8();
});

var count = image_prior.size();
print('Count: ', count);

var palettes = require('users/gena/packages:palettes');
var palette = palettes.matplotlib.magma[7];
var bounds = geometry.bounds()

Map.centerObject(bounds)
var filmArgs = {
    min: 0.00001,
    max: 0.0003,
    maxPixels: 27017280,
    palette : palette,
    opacity:0.75,
    region: geometry,
    crs:'EPSG:3857',
    dimensions: 480,
    framesPerSecond:2,
};

print(image_prior.getVideoThumbURL(filmArgs));


texts: ESRI REST API - Vehicle Routing Problem - Invalid JSON for GPFeatureRecordSetLayer

texts:
I am making a call to the ESRI VRP REST endpoint documented here: http://resources.arcgis.com/en/help/arcgis-rest-api/index.html#/Vehicle_Routing_Problem_service/02r3000000n4000000/

First, despite the url being in the format of a GET request, you have to submit it as a POST request or you get a 403. I have been able to submit a job but my job errors out with the following error:

{
  "type":"esriJobMessageTypeError",
  "description":"Invalid value for parameter Orders - Details : Invalid JSON for GPFeatureRecordSetLayer or GPRecordSet datatype "
}

I have quite a long set of orders, and initially thought that the problem might be that the url is getting cut off, so I moved all the relevant data to url based parameters. In both cases the error is the same. Here is a summary of my order object (the full object can be viewed here: https://pastebin.com/KY6rBT1d) I am including two objects cause they constructed two different ways in my system.

{
    "features": [{
        "geometry": {
            "x": -84.501880199286,
            "y": 38.042953690265
        },
        "attributes": {
            "Name": "C57O1535",
            "ServiceTime": 205,
            "PickupQuantities": "203 203 0 0"
        }
    },
    ...
    {
        "geometry": {
            "x": 4581925.8810913,
            "y": -9402862.099
        },
        "attributes": {
            "Name": "F27",
            "TimeWindowStart1": "08:40:00",
            "TimeWindowEnd1": "08:50:00",
            "MaxViolationTime1": 0,
            "ServiceTime": 6,
            "DeliveryQuantities": "4 0 3 2"
        }
    }
]}

I have checked it's validity using jsonlint in the case of the url, I tried removing the {features: ... } enclosure.


texts: Installed PostGIS extension not listed for database?

texts:
I have a PostgreSQL database with the PostGIS extension enabled.

I know that this extension is enabled, since I am using PostGIS functions, and I can actually see the PostGIS functions.

However, PostGIS is not listed in the installed extensions section, and thus I can't alter it nor uninstall it.

If I try

CREATE EXTENSION postgis;

it tells me that PostGIS is already installed, and that I should uninstall it first. However

DROP EXTENSION postgis CASCADE;

tells me that no PostGIS extension is installed.

What can I do? Any ideas?


texts:
A few months ago I inherited an old (10+ years) Java Desktop application that used maps and simple icons to represent point data to various users. The application was writted in Swing and used Esri MapObjects 2.0 and Geotools 8.3 to receive the map images from a local ASP and a local Arc GIS Map Server. All the maps were Raster. After displaying the maps using some extensions of the Esri Layers point data was represented using icons. Simple, and nothing too fancy.

But the users that now use this application complain about its user-friendliness, having as examples easy to use Google Maps Api and a similar, but online, Web-Map application I created to survey natural disasters and display them in a browser. I suggested we abandoned the Java App, in favor of a web-based one, but this is not possible as clients want to be able to use the application even in a total offline environment (the program gives them the ability to draw features on a featureless map if it cannot connect to the database or the map servers). So I was asked to upgrade it.

I started with Swing, but I came to a dead-end fast. Most of the additional UI elements they request are either abandoned addons or unreliable hacks. So I turned somewhat early to JavaFX and after a few attempts created a basic prototype of the new UI, which was received well.

For that demo, I used the new ARCGis Java SDK, and despite being in beta, it worked really well. Maps, Layers, and a basic image visualization run, and it run well. But the problem was, quite a few of the clients are running on old PCs, a few of them do not even have Direct Draw support, and no OpenGL greater than 1.0. Upgrading is out of the question, so I had to abandon that attempt. I tried using the MapObjects and Geotools 11 in the FX ui, but the SwingNode implementation is terribly unstable. Sometimes the same requests works perfectly, and sometimes the exceptions are piling up. While it is not of this forum to report, the problem lies in EDT and FX interoperability. So much for that in the end.

And so I come to my question. Is anyone aware of a library like the ARCGis SDK or the Geotools library that can work with JavaFX? I need only basic functionality, such as creating a basemap layer from a tile server and the ability to create layers on top of it for simple icons and images. If possible, basic map navigation (pan and zoom) would be welcome, but if not the code to accomplish them is easy to find anyway. I need no complicated like terrain analysis and/or geoprocessing.


texts: How do you extract the data behind a raster histogram

texts:
Bit of a newbie, so this may be a daft question.

I have an .ecw picture - it is has no CRS, as it is a picture taken at ground level with an NIR camera. It was originally converted from a RAW image.

I can load it as a raster image, and generate a histogram. However, I would like to interrogate the data behind the histogram. For example, I might want to take peak height (easy), average, median, max, or min or use a calculation that looks at data spread.

All I can find suggests you use the -c flag in GRASS's r.report (e.g. How to extract pixel value counts from a raster in QGIS?'ve kept the flags the same as this example (apart from the selected image, that is!)

However, when I use it, I keep getting a 'Raster Map Not Found' error. Just trying it with a .tiff map image (which has a CRS), I get the same problem.

Am I missing something blindingly obvious?


texts: OpenLayer's graticule doesn't properly work

texts:
I am making a map of Krasnoyarskiy krai subject of Russia. I have defined a custom projection that should be used for this area using proj4 string. I need to show lat/long grid only for this territory. But in doesn't work as I expected.

At first, even that I set constraints for the projection world extent worldExtent: [72, 40, 120, 85], there are still extra latitude lines (pic #1). Secondly, there are empty parts of longitude lines in the north (pic #1). Finally, if I move my map to the northern areas lines are being built absolutely incorrect. (pic #2).

Here is my HTML:

    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <meta http-equiv="X-UA-Compatible" content="ie=edge">
    <title>graticule</title>
    <link rel="stylesheet" href="./libs/v6.1.1-dist/ol.css">
    <link rel="stylesheet" href="./styles.css">

My JS:

window.onload = init;

function init() {

    proj4.defs('customProj1',"+proj=tmerc +lat_0=0 +lon_0=95 +k=0.998 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs");
    ol.proj.proj4.register(proj4);

    const customProj1 = new ol.proj.Projection({
        code: 'customProj1',
        extent: [-90099540.605703328, -90099540.605703328,
        90099540.605703328, 90099540.605703328],
        worldExtent: [72, 40, 120, 85],
        global: false
    });

    const osmStandartMap = new ol.layer.Tile({
        source: new ol.source.OSM({
            url: "https://{a-c}.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png"
        }),
        visible: true,
        title: 'OSMStandart'
    });

    const graticule = new ol.layer.Graticule({
        strokeStyle: new ol.style.Stroke({
            color: 'rgba(1, 1, 1, 0.9)',
            width: 0.5
        }),
        showLabels: true,
        intervals: [12, 8],
        targetSize: 100,
        wrapX: true,
        visible: true,
        title: 'graticule'

    })

    const map = new ol.Map({
        view: new ol.View({
            center: [25630.78311572224, 7235308.151852912],
            zoom: 8,
            projection: customProj1
        }),
        target: 'js-map'
    });

    const layerGroup = new ol.layer.Group({
        layers: [
            osmStandartMap,
            graticule
        ],
    });
    map.addLayer(layerGroup);
}


texts:
I am suppose to make a map.

The specifications are 1:50k scale centered on a certain point. They want to view a distance of 30 miles from that centered point, but when I scale to 1:50k the map is only about an 8 mile range.

How do I do this in ArcGIS 10?


texts: Loading and rendering shapefile in PyQGIS?

texts:
I wanted to know how qgis reads shapefile. Like i have choose any shapefile from qgis dialog then what are the steps that qgis take to read and render this shapefile data.
Wrather it reads full shapefile data into memory or it opens file connection and reads and renders at runtime?

So far i have found this

layer = QgsVectorLayer(path, name, provider)
if not layer.isValid():
  raise IOError, "Failed to open the layer"

# add layer to the registry
QgsMapLayerRegistry.instance().addMapLayer(layer)

# set extent to the extent of our layer
canvas.setExtent(layer.extent())

# set the map canvas layer set
canvas.setLayerSet([QgsMapCanvasLayer(layer)])

I have opened documentation of each function but i am not successfull in understanding flow of how it is reading shapefile


texts:
Due to performance reasons, I would like to implement some short of tiling + pyramids in a WCS service.

I am not so familiar with WCS, so I was not able to figure out whether it is something supported by the OGC standard.

The use case is to have x/y/t time-series of GeoTiff data (often file sizes are in the scale of GiB), derived from satellite imagery exposed to my client application. I understand that some short of pyramid computation or tiling needs to be done prior to the ingestion in the WCS Server (e.g. GeoServer) for my use case.

Is WCS a feasible approach to implementing the use case described above?